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SUMMARY 

The  goal  of  this  task  was  to  investigate  the  theoretical/mathematical  feasibility  of  using 
femtosecond  (pulsed)  laser-generated  filaments  as  wave  guides  for  radio-frequency  or 
high-power  microwave  transmission.  More  specifically,  we  were  interested  in  the  possible 
use  of  laser-induced  filaments  as  waveguides  in  the  atmosphere. 

We  began  by  first  developing  a  high-resolution  model  for  electromagnetic  wave 
propagation.  This  model  was  initially  derived  from  pre-existing  nonlinear  Schrodinger  (NLS) 
formulations  and  software  left  by  Mary  Potasek,  Sukkeum  Kim,  and  Andrew  Paul.  Our 
reference  equations  extend  the  ’’split-step”  formula  of  G.  P.  Agrawal  with  a  diffraction  term 
used  by  Potasek  as  well  as  Paul.  The  software  was  initially  a  variation  of  FORTRAN 
software  designed  by  Andrew  Paul.  We  improved  that  software  by  increasing  its  resolution 
as  well  as  fixing  various  scaling  errors  that  we  detected.  We  also  supplemented  it  with 
MATLAB  plotting.  We  later  replaced  the  FORTRAN  software  with  MATLAB  equivalents. 
Due  to  a  recurring  problem  with  numerical  instabilities,  this  MATLAB  implementation  was 
ultimately  replaced  by  our  latest  MATLAB  software  which  implements  a  new,  faster 
algorithm. 

With  settings  from  actual  experimental  firings,  we  have  performed  a  series  of  tests 
whereby  various  terms  in  our  expansion  are  either  enabled  or  disabled.  Most  notably,  we 
found  that  with  diffraction  disabled,  our  test  pulse  quickly  collapses  to  a  single  singularity  on 
the  z  axis  due  to  the  non-linear  operator  exploding  at  a  single  time-radius  point,  while  all 
other  values  go  to  near  zero.  In  contrast,  with  diffraction,  the  pulse  tends  to  collapse  more 
gradually  to  a  pair  of  “fins”  off  the  z  axis  before  eventually  collapsing  to  a  pair  of 
singularities  if  the  propagation  is  continued  too  far. 

Through  extensive  testing  of  the  software,  we  determined  that  these  singularity 
conditions  were  primarily  a  result  of  the  fact  that  group-velocity  dispersion  (GVD)  was  not 
included.  So,  in  the  final  months  of  the  task,  we  combined  non-zero  GVD  values  with 
diffraction  and  successfully  modeled  propagation  of  a  Gaussian  pulse  out  to  several  meters 
(propagation  distance)  without  singularities  forming. 

To  see  the  cumulative  evolution  of  the  pulse  over  the  analysis  distance,  the  software 
provides  composite  surface  views.  In  these  views,  we  see  (a)  an  initial  drop,  followed  by 
(b)  a  rise  in  magnitude  at  about  2  to  3  meters,  and  then  (c)  a  gradual  drop  at  larger 
distances.  This  pattern  is  most  pronounced  in  our  “energy  pattern”  depictions  where  we 
model  the  distribution  of  the  total  energy  seen  by  a  target  plane  as  the  pulse  quickly  passes 
through  it.  When  viewed  on  a  target  plane  at  a  given  distance,  the  energy  pattern  appears 
as  a  bright  ring  -  such  that  an  initial  Gaussian  pulse  has  collapsed  to  a  very  thin  cylindrical 
shape. 

The  results  in  this  report  are  numerical  approximations  only  and  actual  experiments  are 
needed  to  validate  our  predictions. 
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1  INTRODUCTION 


We  model  the  evolution  of  a  radially  symmetric  pulse  -  defined  as  a  function  of  radius 
(r),  and  pulse  duration  (t)  -  as  it  propagates  down  a  z  axis. 


Figure  1:  pulse  evolution  model 


Our  formulations  and  software  were  originally  derived  from  a  repository  of  old  computer 
files  -  software,  slide  shows,  and  technical  documents  -  left  by  Mary  Potasek,  Sukkeum 
Kim,  and  Andrew  Paul.  This  includes  a  technical  report  related  to  the  work  of  Kim  and 
Potasek  [3],  Within  this  repository  was  software  developed  by  Andrew  Paul  which  models 
the  propagation  of  laser  pulses  according  to  an  approach  found  in  the  book  by  G.  P. 
Agrawal  [1],  Although  not  directly  solving  the  pertinent  problem,  this  implementation  served 
as  a  good  starting  point  for  the  development  of  our  approach. 

We  also  assembled  published  papers  relating  to  the  mathematical  problem  that  we 
hoped  to  solve.  Most  of  these  papers  do  not  directly  address  the  problem  that  we  are 
solving;  but,  many  contain  fundamental  formulas  and  parameter  settings  that  proved  useful 
in  our  solution.  Many  of  our  parameter  settings  are  derived  from  Sprangle  [8], 

We  used  all  this  material  in  order  to  create  a  reference  formula  guide.  This  guide 
served  as  the  sole  document  from  which  we  created  our  software.  Within  the  guide,  we 
extend  the  “split-step”  formula  of  Agrawal  with  a  diffraction  term  used  by  Potasek  as  well  as 
Paul.  Throughout  all  our  testing,  we  found  that  if  this  diffraction  term  was  disabled,  the 
software  was  generally  an  order  of  magnitude  faster  than  with  diffraction.  Therefore,  for  all 
versions,  there  was  always  a  “fast  model”  (without  diffraction)  as  well  as  a  “slower  model” 
(with  diffraction)  in  the  software. 
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2  SOFTWARE  EVOLUTION  DURING  THE  TASK 

Our  software  evolved  in  3  major  phases:  (A)  upgrade  of  pre-existing  FORTRAN  code 
supplemented  with  MATLAB  plotting,  (B)  complete  conversion  to  MATLAB  (with  options  to 
use  the  MATLAB  Distributed  Computing  Toolbox  for  parallel  processing),  and  (C)  the  latest 
MATLAB  software  which  implements  a  new,  faster  algorithm.  This  development  occurred 
from  December  2005  through  September  2006. 


2.1  Original  FORTRAN  implementations 

Our  software  was  initially  derived  from  the  FORTRAN  “crsplit4”  program  created  by  Andrew 
Paul  based  on  previous  software  from  Mary  Potasek.  This  software  uses  an  iterative  (2- 
part)  symmetric  split-step  approach  with  sparse  coverage.  It  used  FFT  transforms  in  the 
time  direction  only  -  with  the  Numerical  Recipes  FFT.  It  also  used  the  Numerical  Recipes 
tri-diagonal  solver. 


/' 


Dr, 


Dr, 


/ 


m+K) 


Dr, 


Dr, 


Figure  2:  the  iterative  split  step  approach 


Initial  improvements  to  the  software  consisted  of  upgrades  in  source  readability, 
execution  speed,  and  graphical  displays  -  including  the  capability  to  automatically  generate 
MATLAB  display  software  with  the  analysis  FORTRAN  program  at  each  execution. 
Additionally,  parameters  were  modified  to  better  model  propagation  through  air;  Kim  and 
Potasek  had  modeled  the  propagation  of  pulses  in  reverse  saturable  absorbers. 
Subsequently,  various  adjustments  to  the  mathematics  and  associated  software 
implementation  allowed  us  to  improve  the  quality  of  the  answer. 

First,  array  sizes  were  increased  and  z-axis  step  sizes  were  decreased  to  achieve  better 
answers.  This  higher  density  insures  that  our  piece-wise  linear  approximation  of  a  non¬ 
linear  phenomenon  achieves  higher  accuracy.  Also,  in  order  to  minimize  end  effects,  we 
extended  the  time  axis  well  beyond  the  limits  of  interest. 

With  this  denser  coverage,  the  symmetric  split-step  approach  could  safely  be  replaced 
by  a  simple  split-step  approach  that  is  described  in  the  book:  “Nonlinear  Fiber  Optics”  by  G. 
P.  Agrawal.  This  approach  is  faster  and  avoids  dealing  with  any  concern  about 
convergence  of  the  second  (iterative)  part  of  each  symmetric  approach  cycle  --  Paul’s 
software  contained  no  provision  for  insuring  convergence. 
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Figure  3:  the  simple  split  step  approach 


Secondly,  it  was  discovered  that  the  Numerical  Recipes  implementation  of  the  FFT  was 
not  accurate  enough  for  our  purposes.  Simple  tests  involving  repeated  forward  and  reverse 
FFT  cycles  showed  that  the  original  signal  is  not  preserved.  In  these  tests,  the  original 
pulse  was  first  transformed  by  (a)  a  forward  FFT  in  the  time  direction,  and,  then,  (b)  the 
reverse  transform.  This  was  repeated  for  several  passes  -  in  order  to  model  the  effect  of 
the  FFT  on  our  propagation  model.  Therefore,  the  FFT  calls  were  immediately  replaced 
with  calls  to  an  IMSL  implementation  which  does  preserve  signals. 


initial  pulse  after  1  pass  after  10  passes 


Figure  4:  test  result  for  Numerical  Recipes  FFT 

initial  pulse  after  1  pass  after  10  passes 


Figure  5:  test  result  for  IMSL  FFT 


For  quality  assurance,  the  software  has  now  been  converted  to  SI  (MKS)  units  and 
enhanced  with  logic  to  calculate  performance  statistics.  Using  the  new  statistics,  we 
improved  the  speed  of  the  software  as  well  as  its  ability  to  model  propagation  at  finer  z-axis 
increments.  Test  runs  modeling  1500  z-axis  steps  involving  2048  by  512  element  matrices 
have  been  successfully  done  -  requiring  5  hours  on  a  2.8  GHz  computer.  At  that  time, 
shorter  tests  involving  300  steps  required  roughly  1  hour. 

This  software  suffered  from  sensitivities  to  numerical  noise  such  that  it  ended 
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prematurely  for  most  test  runs.  Because  of  problems  with  numerical  instabilities,  the 
software  was  modified  to  use  a  discrete  Hankel  transform  (DHT)  in  the  radial  direction. 
Although  this  was  later  determined  to  be  an  invalid  modification,  it  was  used  in  these  early 
implementations  as  part  of  the  diffraction  calculations. 

Thus,  the  full  implementation  followed  this  sequence  of  calculations  for  each 
propagation  step: 

1.  initial  pulse 

2.  output  of  initial  intensity  and  pulse  energy 

3.  propagation  loop  (per  z  slice) 

(a)  step  to  next  z  (by  constant  increment) 

(b)  forward  FFT  (time) 

(c)  forward  Hankel  Transform  (radius) 

(d)  apply  diffraction  linear  operator  in  (full)  frequency  domain 

Crank-Nicolson  propagation 
tridiagonal  matrix  solver 

(e)  inverse  Hankel  Transform  (radius) 

(f)  apply  dispersion  linear  operator  in  (time)  frequency  domain 

(g)  time  derivatives  of  adjusted  pulse  in  (time)  frequency  domain 

(h)  inverse  FFT  (time) 

(i)  apply  non-linear  operator  in  time  domain 

(j)  output  intensity  and  pulse  energy  (at  specified  z  values) 

In  response  to  questions  about  parallelizability,  we  did  a  thorough  analysis  of  the 
FORTRAN  software  and  concluded  that  the  propagation  phase  (the  time  consuming  part)  of 
the  software  can  be  partially  parallelized  for  significant  performance  improvements.  Timing 
statistics  show  that  3/4  of  the  execution  time  is  spent  in  the  forward/reverse  DHT 
calculations  -  which  are  fully  parallelizable  matrix  operations.  Parts  of  the  linear  and 
nonlinear  operator  applications  are  also  parallelizable. 
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2.2  Early  MATLAB  implementations 

During  April,  we  began  investigating  the  changes  necessary  to  allow  the  code  to  be  used  in 
a  distributed  processing  environment.  Because  of  restrictions  preventing  the  migration  of 
the  IMSL  libraries  to  this  new  system,  as  well  as  the  desire  to  simplify  source  code,  it  was 
decided  to  convert  the  program  completely  to  a  MATLAB  implementation. 

We  successfully  converted  the  existing  FORTRAN  program  into  an  equivalent  MATLAB 
application.  This  new  program  ran  more  quickly  than  the  FORTRAN  version  as  a  result  of 
the  re-expressing  of  time-consuming  parts  as  matrix  operations  (in  place  of  loops)  with  pre¬ 
allocated  arrays.  For  the  same  300-step  calibration  test  case,  the  following  execution  times 
are  seen:  (a)  26  minutes  with  FORTRAN,  versus  (b)  only  4  minutes  with  MATLAB. 

Additionally,  this  new  implementation  requires  fewer  lines  of  code  than  the  FORTRAN 
equivalent  and  takes  advantage  of  the  MATLAB  library  of  specialized  functions  which 
include  Bessel  functions,  FFT  evaluations,  and  trapezoidal  integration.  Therefore,  we 
believe  that  this  translation  benefited  the  quality  of  the  software. 

In  the  process  of  translating  the  software  to  MATLAB,  we  immediately  detected  and 
corrected  several  scaling  errors.  These  errors  involved  the  calculation  of  the  pulse  energy 
plus  actual  and  critical  power.  With  these  fixes,  model  evolution  appeared  to  occur  within  an 
order  of  magnitude  of  matching  the  experimental  observations  2  to  4  meters  versus  20  to 
30  meters. 

We  successfully  added  a  graphical  user  interface  (GUI)  to  the  MATLAB  implementation 
in  order  to  allow  the  user  to  quickly  control  relevant  parameters  for  both  calculations  and 
animations  -  we  enhanced  the  software  by  adding  the  automatic  generation  of  AVI  files 
showing  animated  model  evolutions. 


Figure  6:  MATLAB  Graphical  User  Interface 
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2.3  Attempts  at  Parallelization 

During  May  and  June,  the  software  was  adapted  to  run  with  the  MATLAB  Distributed 
Computing  Toolbox  in  a  limited  parallel-processing  environment.  With  this  system,  matrix 
calculations  are  distributed  from  a  client  machine  to  multiple  “workers”  on  a  “Beowulf” 
cluster  and  processed  in  parallel.  The  new  software  will  also  run  on  a  stand-alone  MATLAB 
version  7  client  (without  connections  to  the  Beowulf)  as  a  conventional  MATLAB 
application. 

Although  our  implementation  is  faster  for  some  test  runs,  the  software  is  actually  slower 
for  the  more  realistic  case  as  data  transfer  penalties  overwhelm  parallelization  savings. 


2.3.1  Success  with  our  Fast  Model  (without  diffraction) 

By  removing  the  diffraction  calculations,  all  array  operations  can  be  sliced  by  groups  of  radii 
such  that  the  propagation  can  be  calculated  with  one  parallel  job. 


•  worker  1  =  (radii  1  to  64)  - 
- propagtion  steps  < - 


>  worker  2  =  (radii  65  to  128) 
- propagtion  steps  < 


□ 


j - ‘worker  3  =  (radii  129  to  192)- 


-  propagtion  steps 


^worker  4  =  (radii  193  to  256)- 
- propagtion  steps  < - 


Figure  7:  Distributed  Processing  without  Diffraction 

With  diffraction  disabled,  the  application  shows  a  modest  speed-up  when  8  workers  on 
the  cluster  are  used  (see  below). 


3.6  GHz  client  only 


3.6  GHz  client  plus 
8  workers  (3.6  GHz) 


500  steps  5000  steps  50000  steps 


client  cpu 

elapsed 

10.1  sec. 

17.9  sec. 

client  cpu 
elapsed 

92.5  sec. 

169.1  sec. 

client  cpu 
elapsed 

820.3  sec. 

1600.6  sec. 

client  cpu 
elapsed 

3.7  sec. 

14.6  sec. 

client  cpu 
elapsed 

1 .2  sec. 

42.7  sec. 

client  cpu 
elapsed 

2.6  sec. 

371.4  sec. 

Figure  8:  Timings  without  Diffraction 
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2.3.2  Problems  with  our  Slower  Model  (with  diffraction) 

As  opposed  to  the  fast  version,  which  uses  a  single  parallel  job,  our  implementation  with 
diffraction  requires  multiple  parallel  jobs  per  propagation  step. 


original  pulse 


J 


(per  propagation  step) 

animation 

~i 

'  —  1  ■ 

1=^1  =  parallel  job  distributed  by  radial  group  (FFT  processing) 

|  =  parallel  job  distributed  by  temporal  group  (DHT  processing) 


Figure  9:  Distributed  Processing  with  Diffraction 

As  a  result,  this  full  version  actually  runs  significantly  slower  (more  than  100  times) 
when  distributed  this  is  apparently  due  to  submission  costs  associated  with  each  parallel 
job  invocation. 


50  steps  500  steps 


client  cpu 

32.3  sec. 

client  cpu 

294.0  sec. 

3.6  GHz  client  only 

elapsed 

33.1  sec. 

elapsed 

300.9  sec. 

3.6  GHz  client  plus 

client  cpu 

58.6  sec. 

8  workers  (3.6  GHz) 

elapsed 

1414.9  sec. 

(est.  4+  hours) 

Figure  10:  Timings  with  Diffraction 

After  contacting  MathWorks  (the  authors  of  MATLAB),  we  received  a  reply  which 
confirmed  that  efforts  to  use  the  MATLAB  version  7  Distributed  Computing  Toolbox  for 
speed  improvements  will  prove  fruitless  as  data  transfer  penalties  (associated  with  the 
transmission  of  megabyte  matrices  from  client  to  worker  set  and  back)  overwhelm 
parallelization  savings;  so,  efforts  to  implement  a  distributed  version  of  the  software  were 
abandoned  -  until  a  new  release  of  their  software  provides  a  satisfactory  option. 
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2.4  Latest  MATLAB  implementation 

In  June,  while  investigating  the  prospect  of  adding  new  terms  (multi-photon  absorption  and 
the  Drude  model),  it  was  determined  that  a  fundamental  error  had  been  made  in  the 
implementation  of  the  diffraction  operator.  Our  approach  to  mix  the  DHT  with  the  Crank- 
Nicolson  propagation  and  its  tri-diagonal  matrix  solver  was  flawed.  When  the  problem  was 
corrected,  the  software  had  difficulty  modeling  more  than  a  few  steps  due  to  numerical 
noise.  Among  other  weaknesses,  the  tri-diagonal  matrix  solver  is  very  sensitive  to 
numerical  errors  due  to  the  fact  that  it  involves  a  lot  of  divisions. 

As  a  result,  a  new,  simpler  approach  was  developed  for  diffraction.  This  led  to  a  major 
revision  of  the  entire  algorithm.  In  the  new  algorithm,  costly  (and  hard  to  debug) 
calculations  involving  exponentiation,  the  Crank-Nicolson  method,  and  a  tri-diagonal  matrix 
solver  were  replaced  with  a  simpler  approach  involving  only  Fourier  and  Hankel 
transformations  to  evaluate  derivatives. 

Although  this  new  algorithm  appeared  to  more  stable  than  earlier  versions,  it  suffered 
from  numerical  noise  problems  at  large  propagation  distances.  So,  the  new  approach  was 
supplemented  to  optionally  use  frequency-weighed  derivative  operators.  These  weightings 
appear  to  filter  out  numerical  noise  such  that  analyses  out  to  5  meters  (involving  20,000 
steps)  can  be  successfully  done  without  any  overflow  problems. 

The  resulting  software  is  very  fast.  We  were  able  to  generate  results  for  each  20,000- 
step  analysis  in  less  than  3  hours  of  software  execution  on  a  single  3.6  GHz  computer.  The 
optional  use  of  the  MATLAB  Distributed  Computing  Toolbox  has  been  included  in  the  new 
software  in  updated  form;  but,  it  is  untested  as  the  software  has  proven  to  be  sufficiently 
fast  on  a  conventional  computer. 

The  following  section  describes  the  mathematics  used  in  our  latest  MATLAB  software. 
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3  LATEST  FORMULA  GUIDE 

In  the  slowly  varying  envelope  approximation,  the  electric  field  is  modeled  as: 


E(z,r,u> )  =  A(r,oj  —  oj0)  exp(+ikoz)  +  A*(r,  oj  +  w0)  exp(— ikoz) 

For  our  analysis,  we  wish  to  model  the  evolution  of  the  field  envelope  (A)  as  the  pulse 
moves  down  a  propagation  axis.  Our  reference  equations  extend  the  “split-step”  formula  of 
Agrawal  with  a  diffraction  term  used  by  Potasek  as  well  as  Paul: 

=  [Dds  +  JV)  A(z,r,t )  +  (A?)  A(z,r,t) 


where 

A  is  the  amplitude  of  the  envelope  of  the  electric  field, 

^ds  is  the  operator  which  accounts  for  absorption  and  dispersion  in  a  linear  medium, 
A  is  the  operator  which  governs  the  nonlinear  effects  induced  by  the  medium,  and 
Ddf  is  the  operator  which  accounts  for  diffraction. 


We  model  the  evolution  of  a  radially  symmetric  pulse  -  defined  as  a  function  of  radius 
(r),  and  pulse  duration  (t)  -  as  it  moves  down  a  propagation  (z)  axis.  The  pulse  is 
normalized  for  numerical  stability  -  we  define: 


Q(z,r,t) 


A(z,  r ,  t) 
.4(0,  0,0) 


A(z,r,t) 

Ao 


The  devisor  is  the  peak  value  of  an  initial  Gaussian  shape. 

The  analysis  appears  sensitive  to  how  we  specify:  (a)  peak  intensity,  (b)  maximum 
power,  (c)  critical  power,  and  (d)  pulse  energy.  Various  published  formulations  are  included 
as  options. 

Also,  we  have  noticed  that  if  diffraction  is  disabled,  the  software  is  generally  an  order  of 
magnitude  faster.  Therefore,  there  is  a  “fast  model”  (without  diffraction)  and  a  “slower 
model”  (with  diffraction). 
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3. 1 1nitial  Gaussian  Pulse 


The  initial  pulse  is  defined  as  a  function  of  scaled  time  (r)  and  scaled  radius  (p). 

t 

Tp 
r 

Tr 

where  tp  is  the  “FWHM”  temporal  width  and  Tr  is  the  “FWHM”  radial  width. 

(Note:  FWHM  denotes  “full-width  half  maximum”). 

An  optional  chirp,  A cp  (applied  to  phase  cp),  is  assumed  to  be  quadratic  in  time,  such  that: 

The  instantaneous  frequency  increases  linearly  from  the 

leading  to  trailing  edge  for  Acp  >  0  which  is  called  “up-chirp” 
while  the  opposite  occurs  for  Acp  <  0  which  is  called  “down- 
chirp”.  If  Acp  =  0,  there  is  no  chirp  and  the  pulse  is 
Gaussian. 

For  our  chirped  Gaussian  pulse,  its  normalized  initial  shape  is: 

T2  T2 

Q(0,r,t)  =  exp(-  —)  exp(-iA^y)  exp 

or  equivalently 

2  /  2 
Q(0,r,t)  =exp[-y(l  +  «A0)]exp  (  -y 


Figure  11:  Gaussian  Pulse 
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3.2  Power  and  Energy 

Power  and  energy  are  calculated  as  follows: 

Let  the  intensity  be: 


I{z,  r,t)  =  \Q(z,r,t)\* 

Then,  total  pulse  energy  [joules]  is: 


(z)  =  J  J  J  r\A(z,r,t)\2d$drdt 
=  27 v  J  j  r\A(z,  r,t)\2drdt 
=  27t^  J  J  rl(z,r,t)drdt 


Figure  12:  cylindrical  integration 

Within  the  FWHM  limits,  the  initial  pulse  energy  for  a  unit  Gaussian  is 

Burzit  =  7T*TR27>(1  -  e_1)er/(  1) 
w  0.9  4  427t(t^)t> 
w  2.9 66(t^)7> 
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Thus,  the  peak  intensity  [Watt/m]  is 


..  a*  __^2_ 


The  peak  intensity  determines  the  maximum  power  [Watt] 


Pmax  =  7T(r|)P0 

which  is  compared  against  the  critical  power  [Watt]  estimate: 
-  Paul’s  formula  (from  FORTRAN  software): 

/  4  U2Uo 

-  formula  from  Sprangle  [8]  : 


■  formula  from  Mechain  [4]: 


Per  it  ( 


Pcrit  —  ( 


2nri2no ' 


3.37A2  , 
8irri2no ' 
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3.3  The  Fast  Model  (without  diffraction) 

Assuming  diffraction  is  negligible, 


d A(z,r,  t ) 
dz 


(pda  +  N)  A(z,  A  t) 


where 

A  is  the  field  envelope, 

Dds  is  a  differential  operator  which  accounts  for  absorption  and  dispersion  in  a  linear 
medium  given  by. 


a  i ,  (2)  <92  1 ,  d3  1 

Dds~~2k  dt2  +  6k  dt3~2a 

and  ^  js  the  nonlinear  operator  which  governs  the  nonlinear  effects  induced  by  the 
medium  -  given  by: 


2irn2TR  d\A 

H - x - ) 


A 


dt 


where: 

k2  is  the  second-order  group  velocity 
dispersion 

(3) 

k  is  the  third-order  group  velocity 
dispersion 

a  is  the  linear  loss  (a<  0  implies  gain), 

D2  is  the  nonlinear  refractive  index, 

A  is  the  central  wavelength  of  the  incident 
laser  pulse, 

c  is  the  speed  of  light  in  the  vacuum, 

Tr  is  the  slope  of  the  Raman  gain  which  is 
related  to  the  delayed  (coherent)  time 
response  of  the  medium.  And 
A  is  the  complex  conjugate  of  A. 
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Equivalently, 


dA 

dz 


i“’SiA 


1 


aA 


ri2  f  .*<9  A  dA* 

-{2Am+AlH 

^2im2TR^d\A\2  A 


A 


dt 


A 


For  a  small  propagation  step  (A z),  approximate: 

dA  AA 
_  _ 

dz  A  z 

A(z  +  A z,r,t)  ~  A(z,r,t)  +  AA 


Then: 


AA 

Az 


~  A  A 

A  dA 

+  A  dA 

1  . 

-  -aA 

_j_  ^  2 7in-2^  |  ^  1 2 


“)  I A  |  A 


^(2A*—  +  A— 

C  ’/  5/:  +  ^  dt 

27rn2TR^d\A\2  ^ 


A 


<9t 


A 
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3.3.1  Spectral  Derivation 

Fourier  transforms  can  be  used  to  evaluate  derivatives: 


Ff  iojtFt 


A 


F 


-1,  ,2 


<4  Ft 


A 


d3 
dt 3 


A  = 


Ff  ico?  Ft 


A 


where 


Ft  denotes  the  Fourier  transform  in  the  temporal  direction,  and 
c°t  is  temporal  frequency. 

With  substitutions: 


AA 


-k{2)  Az 


Ff^Ft 


1 


(-k^Az)  [Ft-HuJ?Ft 


A 

A 


,a 


+  *( 


27rri2Az^ 

A  ; 


IAIM 


|FfWi 


c 

,ri2Az 


A 


)A 


Ff  1iu>tFt 

i-i 


A* 
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Let 


Q(z,r,t) 


A(z,r,t) 
A{  0,0,0) 


A(z,r,t) 

Ao 


Q(z  +  A z,  r,  t)  ~  Q(z,  r,  t)  +  A Q 


A  Q 


AA 

Ao 


Then,  if  Ao  is  real: 

AqAQ  «  i  (^(2)A^) 

-  (i*(3>A*)  [Ft-HojfFt]  A0Q 

-  (|A  z)AoQ 

+  i(^)\A0Q?A0Q 

_  |  (  ’^2^A)Aoq*Aoq  j  [Ft~HujtFt\  A0Q 

-  (i^l)(AoQ)^  [Ff^Ft]  A0Q* 

-  *((2™2JfiAVog)  [Ff'wtFt] \A0Q\2 
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or,  equivalently: 


A Q  k  i  (it(2)Az)  [Fp'wf.F,]  Q 

-  (ifc(3)Az)  [Fr'iwfF,]  Q 

0  L  J 

-  (|a *)g 

-  2((^).43)(Q*)Q[F,-1MFi]« 

-  ((^)A§)  Q2  [Ff  Wt]  Q‘ 

-  %  ( ( TrJ  q  [j?rW,]  \q\2 

where  Q*  is  the  complex  conjugate  of  Q. 
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A  Fast  Fourier  Transform  (FFT)  is  used  here.  Temporal  frequencies  (in  units  of  radians/ 
time)  are  calculated  as: 


A  ujt  = 


27T 

Nt(At) 


71  7T 

-Ut<  At 


The  ordering  of  the  frequencies  follows  standard  FFT  layout. 

wt  =  0  occurs  as  the  first  matrix  element,  and 

negative  frequencies  occur  in  the  second  half  of  the  matrix. 
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3.4  The  Slower  Model  (with  diffraction) 

With  diffraction  included,  a  new  term  is  added  to  the  starting  equation: 


dA, 

dz 

dA 

dz 

dA 

Ih 


\  fast 


fast  + 


fast 


( Ddf )  A(z,r,t) 

ic  (  _  J_  _d\  2 

2no^0  \  ujodt)  ^ r 

iA  ( _ A__<9  \  2 

47rno  \  27 rc  dt )  ^ r 


A 

A 


where 

A  is  the  field  envelope  at  the  beginning  of  the  propagation  step, 
Afast  is  the  field  envelope  calculated  by  the  fast  model, 

Adif  is  the  diffracted  field  envelope, 
no  js  the  linear  refractive  index. 

js  the  central  angular  frequency  of  the  incident  laser  pulse, 
c  is  the  speed  of  light  in  the  vacuum,  and 
A  is  the  central  wavelength  of  the  incident  laser  pulse. 

Integrating 


Adif(z,r,t)  =  Afa3t(z,r,t)  +  J 


iA  / _ \  2 

47T77.0  l  27 rc  dt  J  ^ r 


A 


dz 


As  a  simple  approximation 


Adif(z,r,t )  «  Afast(z,r,t )  + 


iA  f _ A_j9  \  2 

47mo  \  27 rc  dt  /  ^ r 


A 


The  Laplacian  with  azimuthal  symmetry  is  given  by 

2_ld_  dd_ 

r  dr  ^  dr 2 
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So 


A<nf(z,r,t ) 


Afast(z,  r,  t)  +  Az 


iX 

AtTUq 


A  fast  (z  j  T ,  t)  4""  AAdif 


_^d_\  fld_ 
2nc  dt)  V  r  dr 


A 


Where 


AAdif  =  ( 


iXAz ^ 

47mo y 


1  - 


2l—  \  til 

2ttc  dt  J  \  r  dr 


A 


or,  equivalently, 


a  .  AXAz . 

AA^.f  —  (- - ) 

47rnn 


1  d  d2 
r  dr  ^  dr 2 


A  1  d  f  d\  A  d  /  d2 
27 tc  ^r^ dt  \drj  2itdt  \dr2 


A 
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3.4.1  Spectral  Derivation  in  Two  Dimensions 

Fourier  transforms  can  be  used  to  evaluate  derivatives  in  the  radial  and  time  directions: 

—A  =  [f-VvfJ  a 

or  L  J 

^A=-  [F~WrFr]  A 

r\ 

-A  =  [Ft~HujtFt\  A 

where 

Fr  denotes  the  Fourier  transform  in  the  radial  direction  (the  Hankel  transform), 
u r  is  radial  frequency,  and 

Ft  denotes  the  Fourier  transform  in  the  temporal  direction, 
is  temporal  frequency. 

With  substitutions, 

r\  r\ 

-(—A)  =  iFfHuJtFt]  iFAHoJrFr]  A 
ot  or  1  J  L  J 

|(^)  =  -  [*rw<]  [frV2*;]  ^ 

After  rearranging  terms: 

A Adif  =  (- - )(-)  Fr  hujrFr  A 

47 Too  T  L  -l 

-  C-^)[Fr-'w2rFr]A 
47mo  L  J 

+  [™]  ^f‘\a 
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And  the  field  adjustment  for  diffraction  remains: 


Let 


Adif(z  ~  A  fast.  (z ,  r,  t)  T  AAdif 


Q(z,r,t) 


A(z,r,t ) 
4(0, 0,0) 


A(z,r,t) 


Then: 


Qdif(z,r,t)  tv  Qfast(z,r,t )  +  AQdif 


where 


.  ^  .iXAz.,1. 

AQdif  =  (- - )(-) 

47rno  r 

iXAz 


-  (t=t) 


47mo 

M7m0A27rcnr; 
AXAz..  A  .  - 

+  (  —  ^7—) 


Fr  iujrFr 

Q 


Q 


Fr  1U>rFr 


F-  1icurFr 


Ft  1iojtFt  Q 


47TOo  27TC 


F^COrFr 


Ft  licotFt  I  Q 


As  with  the  fast  model,  a  fast  Fourier  transform  (FFT)  is  used  in  the  time  direction  (Ft). 


A  discrete  Hankel  transform  (DHT)  is  used  in  the  radial  direction  (Fr). 
For  this  DHT,  radii  are  non-uniformly  spaced: 

Gj 

r  .  —  F?  - - — 

/  j  ±  l max  ri 

^Nr+ 1 

0  <  j  <  Nr 
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and  radial  frequencies  (in  units  of  radians/distance)  are  also  non-uniformly  spaced: 


(k-V )  j 


Rr. 


-Cj 


0  <  j  <  Nr 

where  Cyare  reference  coefficients  of  the  DHT  implementation. 
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3.4.2  Alternate  weighted-frequency  formulations 

The  central  difference  approximation  to  a  derivative  is: 

d Q  Q(z,r  +  Ar,  t)  —  Q(z,r  —  Ar,  t) 

a7(2’r'()  = - mT - 

This  “leap  frog”  evaluation  can  be  approximated  by  a  ramped  scaling  of  the 
spectral  operator  (in  the  frequency  domain) 

)  FAQ) 


F^-qA)=w 


^ R  \0Jr 

Ur 


where 


Fr  denotes  the  Fourier  transform  in  the  radial  direction  (the  Hankel  transform), 
js  radial  frequency,  and 

u 'R  is  the  maximum  radial  frequency  in  the  transform. 

Let  the  “triangular”  weighting  be  defined  as: 


WR 


Ur  -  |gy 
Ur 


Then,  for  the  “leap  frog”  evaluation: 


dQ 

dr 


[F~1iojWRFr}Q 
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The  diffraction  term  can  be  adjusted  to  include  this  weighting: 


AQd-, 


’if  ~  (!^£)(i)  \F-l^rWRFr]  Q 
47mo  r  1  J 

-  &F^W«F'}Q 
,*AAz)(  A  >(_!) 


Anno  27tc  r ' 


iXAz  A 
^  Anno  '  27 rc 


f; 


Alternately,  more  general  weighting  formulations  may  be  used: 


WR 


Ur  - 


to 


m 

R 


where  m  is  an  integer  exponent  (0,1, 2, 3, 4, 5,...) 


A  “circular”  weighting  is  given  by: 


WR  = 


oX#  —  CO 

UR 


2 

r 


A  “quartic”  weighting  is  achieved  by: 


WR 


Ur  —  OJ 
Ur 


4 

r 


Ft  1icotFt  Q 
1iujtFt\  Q 
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A  “quintic”  weighting  is  achieved  by: 


WR  = 


CO 


R 


CO' 


R 


Again,  the  diffraction  term  is  modeled  as: 


AQ& 


'if 


OA£  1  \F  , 

47m0  r  L  r  J 


AXAz . 

'  A _ _  ' 


F~WrWlFr 


Q 


-  ( 


47mo 
47 mo  27tc  r  L 


AXAz  X 
Attho  27 re 


F-V2W|Fr 


Ft  1iootFt  Q 
1iootFt  \  Q 
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4  PARAMETER  SETTINGS 

Parameters  are  given  to  the  software  through  its  Graphical  User  Interface  (GUI) 
implemented  in  MATLAB  version  7  format. 


4.1  Required  Input  Parameters 

The  following  quantities  (in  MKS  units)  are  specified  at  start-up: 
material  parameters: 

a:  linear  loss  term  (per  meter)  -  note:  negative  implies  gain 

k  :  the  second  order  group  velocity  derivative  (s  /m) 

(note:  negative  implies  anomalous  dispersion; 
while  positive  implies  normal  dispersion) 

(3)  3 

k  :  the  third  order  group  velocity  derivative^  /m) 

no:  the  linear  refractive  index 

r\2'.  the  nonlinear  refractive  index  (m  /Watt) 

Tr:  the  slope  of  the  Raman  gain  (sec) 

Eo:  the  “FWHM”  energy  of  the  initial  laser  pulse  (joule) 

tp\  the  “FWHM”  temporal  width  (sec)  of  the  input  pulse 
tR\  the  “FWHM”  radial  width  (m)  of  the  input  pulse 
A:  the  pulse  wavelength  (m)  for  vacuum 

Acp:  the  chirp  of  the  pulse(>  0  =  up  chirp,  0  =  none,  <  0  =  down  chirp) 

analysis  dimensions: 

1/14:  the  number  of  “FWHM”  pulse  widths  to  define  min.  and  max.  time 

Nt\  the  number  of  time  grid  points 

1/14:  the  number  of  “FWHM”  pulse  widths  to  define  maximum  radius 

Nr.  the  number  of  radial  grid  points 

Zmax  :  the  maximum  propagation  distance  (m) 

Zstep.  the  propagation  step  size  (m) 

(Note:  FWHM  is  the  “full-width  half  maximum”). 
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4.2  Values  Used 

Here  are  our  parameters  based  on  recent  experiments  to  produce  filaments  -  plus 
additional  values  were  taken  from  Sprangle  [8],  Mlejnek  [5],  Moloney[6],  and  Schwarz[7]: 


material  parameters 

a 

m~  1 

0.0 

linear  loss 

k® 

s2/m 

2.2  x  10-29 

second  order  GVD  term 

it® 

s3/m 

0.0 

third  order  GVD  term 

n0 

1.0 

linear  refractive  index 

n2 

m2/W 

1.0  x  10-23 

nonlinear  refractive  index 

Tr 

s 

0.0 

slope  of  Raman  gain 

1 

Deam  parameters 

E0 

joule 

0.015 

initial  pulse  energy 

X 

m 

8.05  x  10“9 

wavelength  of  input  pulse 

Tp 

s 

1.25  x  10“14 

”FWHM”  temporal  pulse  width  (duration) 

tr 

m 

2.5  x  10“3 

”FWHM”  radius  of  pulse 

A  (f) 

0.0 

chirp  of  pulse 

The  following  constant  is  defined  internally: 
Speed  of  light  in  units  of  m/sec: 


c  =  2.99792456  x  108 
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—  time  axis  - 

total  number  of  time  points 

FWHM  temporal  width 

total  number  of  pulse  widths 


-  radial  axis  - 

total  number  of  radial  points 
FWHM  radius 
total  number  of  pulse  widths 

deriv.  weighting: 

- propagation  direction - 

maximim  propagation  distance 

propagation  (z-axis)  step  size 

number  (z-axis)  slice  figures 

-  media  parameters  - 

linear  loss 

2nd  order  GVD 
3rd  order  GVD 
nonlinear  refractive  index 
slope  of  Raman  gain 

—  beam  parameters - 

chirp 

FWHM  pulse  energy 

wavelength 
critical  power: 
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512 


1.25e-014  sec. 
25 

256 

0.0025  meters 
25 

quintic 

5  meters 
0.00025  meters 

50 

0  per  meter 
2.2e-029  secA2An 
0  secA3/m 
1  e-023  mA2Watt 
2e-016  sec. 

0 

0.015  <oule 
8.05e-007  meters 

Sprangle 


e+011  Watt 
e+01 0  Watt 

It  Settings 


5  LATEST  RESULTS 


5.1  Description  of  Output 

The  following  items  are  calculated  at  each  propagation  step: 

intensity  surface  (normalized)  at  each  propagation  step 


I(z,r,t)  =  \Q(z,r,t)\2 

total  pulse  energy  [joules]  at  a  specific  propagation  step 


Etatal  (-2')  27tPo 


rl(z,r,t)drdt 


Integration  is  limited  to  within  the  “FWHM”  limits  and  done  by  the  trapezoidal  method. 

The  MATLAB  software  plots  the  intensity  surface  only  at  intervals  selected  by  the  user; 
we  generally  produced  20  to  30  slides  for  each  20,000  step  analysis. 

To  see  the  cumulative  evolution  of  the  pulse  over  the  analysis  distance,  the  software 
provides  composite  surface  views.  In  these  representations,  the  3-dimensional  pulse 
shapes  must  be  reduced  to  2-dimensional  cross-sectional  representations  where  the  time 
axis  is  eliminated. 

One  obvious  view  is  the  “maxima  silhouette”  produced  by  using  the  pulse  maxima 
across  time  at  each  radius  and  propagation  distance;  this  is  essentially  the  shape  that  a 
target  positioned  down  the  propagation  axis  would  see  coming  at  it. 

An  alternative  view  is  the  time  integral  across  the  pulse  at  each  radius  and  distance. 
This  gives  the  distribution  of  the  total  energy  seen  by  a  target  plane  as  the  pulse  quickly 
passes  through  it.  For  display,  this  surface  is  normalized  relative  to  its  overall  maximum 
value.  We  call  this  our  “energy  pattern”  evolution. 

For  “target  plane”  views,  this  ’’energy  pattern”  at  a  specific  propagation  distance  is 
plotted  in  Cartesian  coordinates  -  because  of  the  inherent  symmetry,  concentric  circular 
shapes  are  seen  in  this  view. 


30 


5.2  Tests  of  the  Fourier  Transformations 

An  initial  concern  was  the  effect  of  repeated  applications  of  Fourier  transformations  (the 
temporal  FFT  and  the  radial  DHT)  on  the  model  -  specifically  where  accumulated 
calculation  errors  would  corrupt  the  answers.  Additionally,  we  needed  to  know  how  the 
alternate  frequency  weighting  would  adjust  the  answers.  So,  we  conducted  a  series  of  tests 
of  the  net  effect  of  repeated  applications  of  the  transforms  (FFT  and  DHT)  with  only 
weighting  included.  For  these  tests,  20,000  propagation  steps  were  modeled;  for  each 
propagation  step,  the  following  sequence  of  transformations  was  calculated  (in  order): 

A)  forward  FFT, 

B)  inverse  FFT, 

C)  forward  DHT, 

D)  apply  weighting  (if  applicable),  and 

E)  inverse  DHT. 


31 


5.2.1  Case  A:  tests  with  no  weighting 

With  no  weighting,  the  transformations  do  not  seem  to  corrupt  the  pulse  noticeably: 


prop,  distance  =  0.0,  pulse  energy  (FWHM)  =  0.014638  J 


prop,  distance  =  5.000  m,  pulse  energy  (FWHM)  =  0.014S38  J 


radius  (mm) 


time 


radius  (mm) 


time  Os 


Figure  14:  Pulse  before  and  after  20,000  transform  steps 
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Figure  15:  maxima  silhouette  evolution  during  20,000  transform  steps 
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integrated  rel.  intensity 


Figure  16:  energy  pattern  evolution  during  20,000  transform  steps 
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5.2.2  Case  B:  tests  with  circular  frequency  weighting 

With  ’’circular”  weighting,  the  pulse  is  quickly  degraded: 


Figure  17:  maxima  silhouette  evolution  during  20,000  transform  steps 


Figure  18:  energy  pattern  evolution  during  20,000  transform  steps 
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5.2.3  Case  C:  tests  with  quartic  frequency  weighting 

With  “quartic”  weighting,  there  is  a  20  percent  loss  of  signal  during  our  20,000  step 
propagation  test. 


Figure  19:  maxima  silhouette  evolution  during  20,000  transform  steps 


Figure  20:  energy  pattern  evolution  during  20,000  transform  steps 
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5.2.4  Case  D:  tests  with  quintic  frequency  weighting 

With  “quintic”  weighting,  there  is  less  than3percent  loss  of  signal  during  our  20,000  step 
propagation  test. 


Figure  21:  maxima  silhouette  evolution  during  20,000  transform  steps 
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Figure  22:  energy  pattern  evolution  during  20,000  transform  steps 
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5.3  Propagation  Tests  without  Diffraction 

Without  diffraction,  the  pulse  tends  to  collapse  along  the  propagation  axis. 

5.3.1  Nonlinear  terms  only 

With  only  the  nonlinear  terms  active  (the  linear  operator  set  to  zero),  the  pulse  quickly 
evolves  to  a  singular  spike  on  the  propagation  axis. 


Figure  23:  pulse  at  2  and  2.5  meters 


maxima  of  propagated  pulse 


Figure  24:  maxima  silhouette  up  to  the  singularity 


37 


distance  (m) 
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normalized  (tune)  integrated  propagated  pulse 


normalized  (time)  integrated  propagated  pulse 


distance  (m) 


Figure  25:  energy  pattern  evolution  up  to  the  singularity 


5.3.2  Nonlinear  terms  plus  second-order  GVD 

By  contrast,  with  the  second  order  group  velocity  dispersion  (GVD)  non-zero,  the  pulse  is 
stretched  out  in  the  propagation  direction. 


Figure  26:  pulse  shapes  at  0.0, 1.5,  3.0,  and  4.5  meters 


Figure  27:  maxima  silhouette  evolution 
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normalized  (time)  integrated  propagated  pulse 


normalized  (time)  integrated  propagated  pulse 


5.4  Propagation  Tests  with  Diffraction  but  No  Linear  Operator 

We  performed  tests  where  the  linear  dispersion  operator  was  set  to  zero;  but,  diffraction 
was  enabled.  These  tests  combined  a  non-zero  nonlinear  operator  with  the  active 
diffraction. 


5.4.1  Nonlinear  only  with  diffraction  and  no  weighting 

At  a  very  short  propagation  distance,  the  pulse  evolves  to  a  pair  of  singular  spikes 
alongside  the  propagation  axis;  on  the  target  plane,  this  will  actually  be  a  ring. 


Figure  29:  pulse  shapes  at  0.0,  0.5, 1.0,  and  1.5  meters 


Because  the  pulse  goes  singular  at  1 .5  meters,  the  analysis  can  not  continue  to  larger 
propagation  distances. 
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5.4.2  Nonlinear  only  with  diffraction  and  quintic  frequency  weighting 

The  incorporation  of  the  alternate  frequency  weighting  did  not  provide  much  significant 
change  to  our  results. 


prop,  distance  -  600.000  mm.  puke  energy  (FWHM)  -  0.014327  J 


radius  (mm)  fa) 


prop,  distance  -  1.600  m.  pulse  energy  (FWHMt  -  0.013236  J 
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radius  (mm)  time  (fs) 


prop,  distance  -  1.000  m.  pulse  energy  (f  WHM)  -  0.013044  J 


radius  (mm)  time  (Is) 

prop,  distance  -  1.700  m,  pulse  enorgy  (f  WHM)  -  04)12691  J 


"I 

15  J 


radius  (mm)  *imB  (M 


Figure  30:  pulse  shapes  at  0.5, 1.0, 1.5,  and  1.7  meters 

Analyses  of  the  results  of  these  tests  led  us  to  conclude  that  a  non-zero  second-order 
GVD  is  required  if  the  pulse  evolution  is  to  be  modeled  at  larger  propagation  distances. 
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5.5  Propagation  with  Diffraction  and  Non-Zero  Second-Order  GVD 

Our  best  results  appear  to  come  from  models  where  diffraction  and  the  second-order 
GVD  are  both  non-zero  and  frequency  weighting  is  used. 

5.5.1  Singularity  Problems  with  Unweighted  Diffraction 

Without  any  frequency  weighting,  the  model  progresses  up  to  2  meters  with  out  any 
numerical  problems.  However,  at  around  3  meters,  the  model  becomes  corrupted  by 
numerical  noise. 


Figure  31:  pulse  shapes  at  0.5, 1.0, 1.5,  and  2.0  meters 
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prop,  distance  =  2.500  m,  pulse  energy  (FWHM)  =  0.006792  J 


prop,  distance  -  2.900  m.  puke  energy  (FWHM)  -  0.005620  J 


prop,  distance  -  3.000  m.  pulse  energy  (FWHM)  -  0.005354  J 


prop,  distance  -  3.100  m.  puke  energy  (FWHM)  -  0.005100  J 


(overflow  at  3.1975  m.) 


Figure  33:  pulse  corruption  around  3  meters 


Figure  32:  pulse  shape  at  2.5  meters 


prop,  distance  =  2.500  m,  pulse  energy  (FWHM)  =  0.006792  J 


5  0-5 

radius  (mm) 


5.5.2  Best  Results  with  Quartic  Weighting 

By  using  the  “quartic”  frequency  weighting,  we  can  model  the  pulse  beyond  3  meters.  For 
our  tests,  we  managed  to  model  the  evolution  out  to  5  meters. 


prop,  distance  ■  2.000  m,  pulse  energy  (FWHM)  -  0 .000486  J 


prop,  distance  -  2.001)  m.  pulse  energy  (FWHM)  -  0.008480  J 


Figure  34:  pulse  shapes  at  2  and  2.5  meters 
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Figure  35:  pulse  shapes  at  3  and  3.5  meters 


Figure  36:  maxima  silhouette  evolution 
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Figure  37:  energy  pattern  evolution 


prop,  distance  =  2.500  m,  pulse  energy  i'FWHM)  =  0.006857  J 
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Figure  38:  energy  pattern  as  seen  by  target  plane  at  2.5  meters 
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5.5.3  Best  Results  with  Quintic  Weighting 

As  a  result  of  the  findings  of  the  transform  tests,  we  replaced  all  the  “quartic-weighted” 
analyses  with  “quintic-weighted”  equivalents.  This  new  weighting  improves  the  quality  of  the 
answer  beyond  3  meters.  As  with  the  “quartic-weighted”  tests,  we  managed  to  model  the 
evolution  out  to  5  meters. 


piup.  ilhUaiiur  -  ?.000  in.  pub*  ntnriijy  (FW1IM)  -  01004 70  J 


piup.  t  Istantt*  -  2.000  m.  puke  eneigy  (TWHM)  -  0000426  J 


radius  (mm)  W 


piop.  distance  -  2.500  m.  puke  energy  (TWHM)  -  0006831  J 


radius  (mm)  0mo  (fs) 


radius  (mm) 


Figure  39:  pulse  shapes  at  2  and  2.5  meters 
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maximum  rel.  intensity 


Figure  40:  pulse  shapes  at  3  and  3.5  meters 


Figure  41:  maxima  silhouette  evolution 


49 


3  04 

* 

f  0.2 

E 


0 

•10 

radius  (mm) 


distance  (m) 


radius  (mm) 


4 

distance  (m)  5 


Figure  42:  energy  pattern  evolution 


Figure  43:  energy  pattern  as  seen  by  target  plane  at  2.7  meters 


50 


Figure  44:  comparison  of  composite  pulse  evolution  profiles 
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5.5.4  Tests  with  Higher  GVD 

If  the  second-order  GVD  value  is  doubled,  the  following  results  are  obtained: 


Figure  45:  maxima  silhouette  evolution 


Figure  46:  energy  pattern  evolution 
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radius  (mm)  maximum  rel.  intensity 


Figure  47:  comparison  of  composite  pulse  evolution  profiles 


Figure  48:  energy  pattern  as  seen  by  target  plane  at  3.6  meters 
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5.6  The  Effect  of  the  Chirp  Parameter 

An  optional  chirp,  A cp  (applied  to  phase  (p),  is  assumed  to  be  quadratic  in  time.  By  our  convention, 
based  on  that  of  Agrawal,  the  instantaneous  frequency  increases  linearly  from  the  leading  to  trailing 
edge  for  A  cp>  0  which  is  called  ”up-chirp”;  while  the  opposite  occurs  for  Acp  <  0  which  is  called 
”down-chirp”.  If  A  cp  =  0,  there  is  no  chirp  and  the  pulse  is  Gaussian. 

We  repeated  our  “quintic- weighted”  analyses  with  various  values  of  non-zero  chirp. 


5.6.1  Example  with  Positive  Chirp 

A  simple  test  of  positive  chirp,  where  Acp  =  +1.0,  produced  the  following: 


radius  (mm) 


time  (fs) 


prop,  distance  =  2.500  m,  pulse  energy  (FWHM)  =  0.006831  J 


prop,  distance  =  2.500  m,  pulse  energy  (FWHM)  =  0.006170  J 


with 


without  chirp 


radius  (mm) 


time  (fs) 


Figure  49:  pulse  at  2.5  m.  (a)  without  and  (b)  with  chirp  of  +1.0 


radius  (mm) 


time  (fs) 


radius  (mm) 


time  (fs) 


prop,  distance  =  3.000  m,  pulse  energy  (FWHM)  =  0.005414  J 


prop,  distance  =  3.000  m,  pulse  energy  (FWHM)  =  0.004901  J 


without 


Figure  50:  pulse  at  3.0  m.  (a)  without  and  (b)  with  chirp  of  +1.0 
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maximum  rel.  intensity  maximum  rel.  intensity 


Figure  51:  maxima  silhouette  (a)  without  and  (b)  with  chirp  of  +1.0 


maxima  of  propagated  pulse  maxima  of  propagated  pulse 


distance  (m)  distance  (m) 


Figure  52:  maxima  silhouette  (a)  without  and  (b)  with  chirp  of  +1.0 
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Figure  53:  energy  pattern  (a)  without  and  (b)  with  chirp  of  +1.0 


distance  (m) 


distance  (m) 


Figure  54:  energy  pattern  (a)  without  and  (b)  with  chirp  of  +1.0 
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The  peak  in  the  energy  pattern  moves  from  2.7  to  2.9  meters  when  a  chirp  of  +1.0  is  used;  but  a 
similar  pattern  is  seen. 


Figure  55:  peak  energy  pattern  as  seen  by  target  (a)  without  and  (b)  with  chirp  of  +1.0 
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5.6.2  Example  with  Negative  Chirp 

A  simple  test  of  negative  chirp,  where  A <p  =  1.0,  produced  the  following: 


prop,  distance  =  2.500  m,  pulse  energy  (FWHM)  =  0.006831  J 


without  chirp 


radius  (mm)  t'me  (k) 


prop,  distance  =  2.500  m,  pulse  energy  (FWHM)  =  0.007672  J 

with  chirp  of -1  !P"T  T  ^ 


radius  (mm)  time  (fs) 


Figure  56:  pulse  at  2.5  m.  (a)  without  and  (b)  with  chirp  of  -1.0 


prop,  distance  =  3.000  m,  pulse  energy  (FWHM)  =  0.005414  J 


without  chirp 


prop,  distance  =  3.000  m,  pulse  energy  (FWHM)  =  0.006108  J 


time 


radius  (mm) 


time  (fs) 


radius  (mm) 


Figure  57:  pulse  at  3.0  m.  (a)  without  and  (b)  with  chirp  of  -1.0 
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maximum  rel.  intensity  maximum  rel.  intensity 


maxima  of  propagated  pulse  maxima  of  propagated  pulse 


without  chirp  ..  ”  :  with  chirp  of -I 


Figure  58:  maxima  silhouette  (a)  without  and  (b)  with  chirp  of  -1.0 


maxima  of  propagated  pulse  maxima  of  propagated  pulse 


distance  (m)  distance  (m) 


Figure  59:  maxima  silhouette  (a)  without  and  (b)  with  chirp  of  -1.0 
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normalized  (time)  integrated  propagated  pulse 


normalized  (time)  integrated  propagated  pulse 


Figure  60:  energy  pattern  (a)  without  and  (b)  with  chirp  of  -1.0 


Figure  61:  energy  pattern  (a)  without  and  (b)  with  chirp  of  -1.0 
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The  peak  in  the  energy  pattern  moves  from  2.7  to  2.6  meters  when  a  chirp  of  -1.0  is  used;  but  a 
similar  pattern  is  seen. 


radius  (mm)  5  -5 


radius  (mm) 


radius  (mm)  5  -5 


radius  (mm) 


prop,  distance  =  2.700  m,  pulse  energy  (FWHM)  =  0.006231  J 


without  chirp 


prop,  distance  =  2.600  rn,  pulse  energy  (FWHM)  =  0.007332  J 


with  chirp  of-1 


Figure  62:  peak  energy  pattern  as  seen  by  target  (a)  without  and  (b)  with  chirp  of  -1.0 
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maximum  rel.  intensity  maximum  rel.  intensity 


5.6.3  Example  with  Extreme  Negative  Chirp 

A  more  extreme  test  of  negative  chirp,  where  A <p  =  2.0,  produced  the  following: 


maxima  of  propagated  pulse 


maxima  of  propagated  pulse 


with  chirp  of  -2 


radius  (mm) 


distance  (m) 


Figure  63:  maxima  silhouette  (a)  without  and  (b)  with  chirp  of  -2.0 


maxima  of  propagated  pulse 


without  chirp 


1.5  2  2.5  3 

distance  (m) 


maxima  of  propagated  pulse 


Figure  64:  maxima  silhouette  (a)  without  and  (b)  with  chirp  of  -2.0 
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normalized  (time)  integrated  propagated  pulse 


normalized  (time)  integrated  propagated  pulse 


Figure  65:  energy  pattern  (a)  without  and  (b)  with  chirp  of  -2.0 


Figure  66:  energy  pattern  (a)  without  and  (b)  with  chirp  of  -2.0 
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The  peak  in  the  energy  pattern  moves  from  2.7  to  2.4  meters  when  a  chirp  of  -2.0  is  used;  but  a 
similar  pattern  is  seen. 


prop,  distance  =  2.700  m,  pulse  energy  (FWHM)  =  0.006231  J 
without  chirp  ..  B 


& 


prop,  distance  =  2.400  m,  pulse  energy  (FWHM)  =  0.009119  J 
with  chirp  of -2 


* 


Figure  67:  peak  energy  pattern  as  seen  by  target  (a)  without  and  (b)  with  chirp  of  -2.0 
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CONCLUSIONS 

In  our  latest  views,  we  see  (a)  an  initial  drop,  followed  by  (b)  a  rise  magnitude  at  about  2 
to  3  meters,  and  then  (c)  a  gradual  drop  at  larger  distances.  This  pattern  is  most 
pronounced  in  our  “energy  pattern”  depictions  where  we  model  the  distribution  of  the  total 
energy  seen  by  a  target  plane  as  the  pulse  quickly  passes  through  it.  When  viewed  on  a 
target  plane  at  a  given  distance,  the  energy  pattern  appears  as  a  bright  ring  -  such  that  an 
initial  Gaussian  pulse  has  collapsed  to  a  very  thin  cylindrical  shape. 

Our  results  are  based  solely  on  mathematical  formulations  without  any  experimental 
verification.  In  the  future,  we  hope  that  observations  of  actual  laser  experiments  will 
provide  results  that  match  our  predictions. 

Additionally,  the  software  currently  lacks  adequate  safeguards  to  completely  ensure 
energy  conservation  with  the  propagated  pulse.  Future  upgrades  should  include 
diagnostics  and/or  corrections  to  address  this  shortcoming. 
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Appendix  A:  HOW  TO  USE  THE  LATEST  MATLAB  VERSION 

The  latest  implementation  requires  MATLAB  version  7  to  be  installed  on  the  computer  where  the 
application  is  to  be  used.  This  MATLAB  application  is  started  with  a  special  DOS  batch  file: 


[^3  nu  d  Iig  1  iRe  s  NL  S  Jt>  at 
Figure  68:  starting  the  application 

This  causes  the  initial  pop-up  to  appear: 


Figure  69:  the  graphical  user  interface  (GUI)  at  startup 
Calculation  is  started  with  the  “Calculate”  button  at  bottom  center. 
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When  calculations  are  activated,  the  buttons  at  the  bottom  turn  red;  a  step  counter  also  appears  at  the 
bottom  left  to  indicate  status. 


Figure  70:  the  graphical  user  interface  (GUI)  during  active  calculations 
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The  button  turns  back  to  green  when  calculations  are  completed;  plot  buttons  activate  various  display 
capabilities. 


Figure  71:  the  graphical  user  interface  (GUI)  at  completion 
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rel.  intensity 


The  first  plot  button  activate  pulse  animation  displays 


1  Figure  1 

File  Edit  View 

Insert  Tools  Desktop 

Window  Help 

D£B3 

Q  QO  $ 

□  Id 

“bO 

prop,  distance  =  0.0,  pulse  energy  (FWHM)  =  0.014638  J 


1 

0.3 

0.3 

0.4 

0.2 

1 


Figure  72:  pulse  animation  plotting 
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animation  option  C  uncompressed  C  compressed  (•  off 
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close  all  figures 
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The  second  plot  button  activate  composite  profile  displays 


)  Figure  1 

mxi 

File  Edit  View 

Insert  Tools  Desktop 

Window  Help 

'SI 

v-©| 

□  0 

□  □ 

maxima  of  propagated  pulse 


Tb-y  " 

* 


radius  (mm)  distance  (m) 


Figure  73:  composite  profile  plotting 
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norm,  integrated  intensity 


The  third  plot  button  activate  target-plane  energy  pattern  displays 


•}  Figure  26 


File  Edit  View  Insert  Tools  Desktop  Window  Help 

Di£y#|  k|Q.Q.3n)$|sl!|[D'E|  □  0 


prop,  distance  =  2.500  m,  pulse  energy  (FWHM)  =  0.008720  J 
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frame  51  of  51 


close  all  figu-es 


Figure  74:  display  of  energy  pattern  on  target  planes 
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Appendix  B:  LATEST  MATLAB  SOURCE  CODE 


The  MATLAB  software  used  in  this  project  consists  of  two  parts:  (A)  figure  files  for  the  graphical 
user  interface,  and  (B)  executable  MATLAB  textural  source  code. 

B.1  Figures  for  the  Graphical  Interface 

The  graphical  user  interfaces  for  this  program  were  created  with  GUIDE,  the  MATLAB  graphical 
user  interface  development  environment.  GUIDE  creates  two  files,  a  FIG-file  and  an  M-file.  The 
FIG-file,  with  extension  .fig,  is  a  binary  file  that  contains  a  description  of  the  layout.  The  M-file, 
with  extension  .m,  contains  the  code  that  controls  the  GUI  and  responds  to  user  actions. 

B.1.1  Interface  for  the  main  program  (highResNLS) 

First,  there  is  the  interface  for  the  main  program,  highResNLS.fig: 


Figure  75:  GUI  definition  for  the  main  program 
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The  following  identifiers  and  callback  definitions  must  be  defined: 


string 

tag 

style 

callback(s) 

512 

timeGrid512 

radiobutton 

timeGrid512  Callback 

1024 

timeGridl024 

radiobutton 

timeGridl024  Callback 

1.25e-014 

taup 

edit 

taupCallback,  taupCreateFcn 

25 

widths 

edit 

widthsCallback,  widthsCreateFcn 

256 

radialGrid256 

radiobutton 

radialGrid256  Callback 

512 

radialGrid512 

radiobutton 

radialGrid512  Callback 

0.0025 

tauR 

edit 

taurCallback,  taurCreateFcn 

25 

RWidth 

edit 

RWidthCallback,  RWidthCreateFcn 

none 

derivFlagRawR 

radiobutton 

derivFlagRawR  Callback 

quintic 

derivFlagQuinticR 

radiobutton 

derivFlagQuinticR  Callback 

quartic 

derivFlagQuarticR 

radiobutton 

derivFlagQuarticR  Callback 

5 

maxPropDist 

edit 

maxPropDistCallback,  maxPropDistCreateF cn 

0.00025 

dz 

edit 

dzCallback,  dzCreateFcn 

50 

NZSlice 

edit 

NZSliceCallback,  NZSliceCreateFcn 

0 

linearLoss 

edit 

linearLossCallback,  linearLossCreateFcn 

2nd  order  GVD 

gvd2onoff 

checkbox 

gvd2onoffCallback 

2.1e-030 

gvd2deriv 

edit 

gvd2derivCallback,  gvd2derivCreateFcn 

0 

gvd3deriv 

edit 

gvd3  derivCallback,  gvd3  derivCreateFcn 

1.0 

linearRefract 

edit 

linearRefract  Callback, 

linearRefract  CreateFcn 

le-023 

nonlinRefract 

edit 

nonlinRefractCallback,  nonlinRefractCreateF cn 

slope  of . . . 

ramanOnOff 

checkbox 

ramanOnOffCallback 

2e-016 

ramanSlope 

edit 

ramanSlopeCallback,  ramanSlopeCreateFcn 

diffraction 

diffraction 

checkbox 

diffraction  Callback 

0 

chirp 

edit 

chirpCallback,  chirpCreateFcn 

1.5e-002 

pulseEnergy 

edit 

pulseEnergyCallback,  pulseEnergyCreateFcn 

8.05e-007 

wavelength 

edit 

wavelengthCallback,  wavelengthCreateFcn 

Paul 

powerPaul 

radiobutton 

powerPaulCallback 

Sprangle 

powerSprangle 

radiobutton 

powerSprangleCallback 

Mechain 

powerMechain 

radiobutton 

powerMechainCallback 

zoomed? 

zoomOption 

checkbox 

zoomOption  Callback 

distributed? 

parallel 

checkbox 

parallel  Callback 

Calculate 

calculate 

pushbutton 

calculateCallback 
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Show  Animation 

plotl 

pushbutton 

plotl  Callback 

Show 

Cummulative 

plot2 

pushbutton 

plot2  Callback 

Show  Target  Plane 

plot3 

pushbutton 

plot3  Callback 

status 

text 

(none) 

maximum  power... 

maxPower 

text 

(none) 

critical  power... 

criticalPower 

text 

(none) 
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B.1.2  Interface  for  pulse  animation  plotting  (highResNLSplotl) 

For  pulse  animation  plotting,  the  software  uses  highResNLSplotl.fig: 


Figure  76:  GUI  definition  for  the  pulse  animation  plotter 

The  following  identifiers  and  callback  definitions  must  be  defined: 


string 

tag 

style 

callback(s) 

incoming 

incoming 

radiobutton 

incoming  Callback 

outgoing 

outgoing 

radiobutton 

outgoingCallback 

target 

target 

radiobutton 

targetCallback 

uncompressed 

uncompressed 

radiobutton 

uncompressedCallback 

compressed 

compressed 

radiobutton 

compressedCallback 

off 

noanimation 

radiobutton 

noanimation  Callback 

half  view 

halfview 

checkbox 

halfview  Callback 

jpeg  export 

export 

checkbox 

exportCallback 

plot 

plot 

pushbutton 

plot  Callback 

close  all  figures 

closeall 

pushbutton 

closeall  Callback 

status 

text 

(none) 
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B.1.3  Interface  for  composite  profile  plotting  (highResNLSplot2) 

For  the  composite  profile  displays,  the  software  uses  highResNLSplot2.fig: 


Figure  77:  GUI  definition  for  the  composite  profile  plotter 
The  following  identifiers  and  callback  definitions  must  be  defined: 


string 

tag 

style 

callback(s) 

incoming 

incoming 

radiobutton 

incoming  Callback 

outgoing 

outgoing 

radiobutton 

outgoingCallback 

target 

target 

radiobutton 

targetCallback 

maximum 

maxVal 

radiobutton 

maxVal  Callback 

integrated 

integrated 

radiobutton 

integrated  Callback 
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B.1.4  Interface  for  target-plane  animation  plotting  (highResNLSplot3) 

For  animations  involving  target  plane  energy  patterns,  the  software  requires  highResNLSplot3.fig: 


Figure  78:  GUI  definition  for  the  target-plane  animation  plotter 

The  following  identifiers  and  callback  definitions  must  be  defined: 


string 

tag 

style 

callback(s) 

maximum 

maxVal 

radiobutton 

maxVal  Callback 

integrated 

integrated 

radiobutton  integrated  Callback 

uncompressed 

uncompressed 

radiobutton 

uncompressedCallback 

compressed 

compressed 

radiobutton 

compressedCallback 

off 

noanimation 

radiobutton 

noanimation  Callback 

2-D 

draw2D 

radiobutton 

draw2D  Callback 

3-D 

draw3D 

radiobutton 

draw3D  Callback 

plot 

plot 

pushbutton 

plot  Callback 

close  all  figures 

closeall 

pushbutton 

closeall  Callback 

status 

text 

(none) 
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B.2  Textural  Source  Code 


The  figure  files  are  combined  with  MATLAB  textural  source  code  in  order  to  define  the  complet 
executable  software. 

B.2.1  The  main  program  (highResNLS) 


First,  there  is  the  main  program,  highResNLS.m: 

function  varargout  =  highResNLS (  varargin  ) 

9 ''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 
o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  client-side  main  program  with  GUI 

O, 

o 

%  14  Sep  2006  version  5.5 

%  with  option  for 

%  MATLAB  7  Distributed  Computing 

O, 

o 

%  This  program  is  a  simplification  of  a  rewrite  of  crsplit4  by  A.  Paul. 

%  which  is  partially  based  on  the  Split  Step  method  described  in  the 
%  book: 

%  Nonlinear  Fiber  Optics  by  G.P.  Agrawal . 

O, 

o 

%  This  method  for  the  extended  nonlinear  Schrodinger  equation  (NLS) 

%  models  the  propagation  of  optical  pulses  in  a  nonlinear  material. 


input  files 


c  .mat 


--  coefficients  for  DHT 


output  files 


plott .mat 
plotr . mat 
distances . mat 
energies . mat 
intensity . mat 


time  axis  grid 
radial  axis  grid 
propagation  distances 
pulse  energies 
intensity  surfaces 


Control  Parameters  (inputs) 


media  parameters: 

linear  loss  term  per  meter  [1/m] 
2nd-order  group  vel.  deriv.  [secA2/m] 
3rd-order  group  vel.  deriv.  [secA3/m] 
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o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 

g, 

o 

g, 

o 

o, 

o 

g, 

o 

o, 

o 


linear  refractive  index  (unitless) 
nonlinear  refractive  index  [mA2/Watt] 
slope  of  the  Raman  gain  [sec] 
diffraction  on/off  flag 
beam  parameters : 

chirp  (-1  =  down  chirp,  +1  =  up  chirp,  0  =  No  chirp) 
pulse  energy  [joule] 
wavelength  for  vacuum  [m] 
self-f ocusing  power  option 

1  =  Paul's  formula  (AFRL  Tech.  Report) 

2  =  Sprangle  formula  (Physical  Review  E  66,  046418  [2002]) 

3  =  Mechain  formula  (Teramobile) 
time-axis  parameters: 

pulse  width  [sec] 
total  number  of  time  points 
total  number  of  pulse  widths 
radial-axis  parameters: 

FWHM  radius  [m] 
total  number  of  radial  points 
total  number  of  radial  pulse  widths 
derivative  option 
propagation-axis  parameters: 

maximim  propagation  distance  [m] 
propagation  (z-axis)  step  size  [m] 
number  (z-axis)  slice  figures 
output : 

zoom  option: 

0  for  wide  view  (512x256  plotted  with  4x4  averaging) 

1  for  zoomed  view  (center  128x64  plotted) 


Here  are  a  few  parameters  from  recent  experiments  to  produce  filaments: 

Wavelength:  805  nm  (Center)  +/-  10  nm  Bandwidth 

Beam  Radius:  ~15  mm 

Pulse  Energy:  430mJ 

Pulse  Duration:  50fs 

Chirp:  No  chirp 

Beam  Collimated  (non-f ocused) 

Filament  in  2-3  meters  from  output  aperture 
Number  of  filaments  (60  +)  (AKA  A  BUNCH) 

Difficult  to  tell  if  pattern  remains  the  same  shot-to-shot 


Wavelength:  805  nm  (Center)  +/-  10  nm  Bandwidth 

Beam  Radius:  ~5  mm 

Pulse  Energy:  15mJ  (8  mj  is  threshold  for  occasional  filament) 

Pulse  Duration:  50fs 

Chirp:  No  Chirp 

Beam  Collimated  (non-f ocused) 

Filament  in  4  meters  from  output  aperture 
Number  of  filaments:  1 


The  original  software  has  been  extended  to  use  the  Hankel  Transform  of 
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%  M.  Guizar-Sicairos  and  J.  C.  Gutierrez-Vega  --  which  implements  Hankel 
%  transforms  of  integer  order  based  on  a  Fourier-Bessel  series  expansion 
%  as  described  in  the  recently  published  work: 

%  M.  Guizar-Sicairos  and  J.  C.  Gutierrez-Vega,  Computation  of 

%  quasi-discrete  Hankel  transforms  of  integer  order  for  propagating 

%  optical  wave  fields,  J.  Opt.  Soc.  Am.  A  21,  53-58  (2004) . 

%  The  numerical  method  features  great  accuracy  and  is  energy  preserving  by 
%  construction,  it  is  especially  suitable  for  iterative  transformation 
%  processes.  Its  implementation,  requires  the  computation  of  zeros  of 
%  the  m-th  order  Bessel  function  of  the  first  kind  where  m  is  the 
%  transformation  order.  An  array  of  the  first  3001  Bessel  functions  of 
%  order  from  zero  to  four  can  be  found  in  the  "c.mat"  array.  If  a  greater 
%  transformation  order  is  required  the  zeros  may  be  found  numerically. 

%  With  the  c.mat  array,  as  included,  Hankel  transforms  of  order  0-4  may  be 
%  computed,  with  up  to  3000  sampling  points. 

O, 

o 

o _ 

O - 

%  main  program 

O _ 

O - 


format  long; 

%  Begin  initialization  code  -  DO  NOT  EDIT 

gui  Singleton  =  1; 

gui_State  =  struct (  'gui_Name', 

'gui  Singleton', 

'gui  OpeningFcn ' , 

' gui_OutputFcn ' , 

' gui_LayoutFcn ' , 

' gui_Callback ' , 

if (  nargin  &  isstr (varargin { 1 } )  ) 

gui  State. gui  Callback  =  str2func(  varargin{l}  ); 

end 

if (  nargout  ) 

[varargout { 1 : nargout } ]  =  gui  mainfcn (  gui  State,  varargin) : }  ); 

else 

gui  mainfcn (  gui  State,  varargin) : }  ); 

end 

%  End  initialization  code  -  DO  NOT  EDIT 


mf ilename, 
gui  Singleton, 
@highResNLS  OpeningFcn, 
@highResNLS  OutputFcn, 
[], 

[] 


executes  just  before  highResNLS  is  made  visible. 


function  highResNLS_OpeningFcn (  hObject,  eventdata,  handles,  varargin  ) 


%  hObject 
%  eventdata 
%  handles 


handle  to  figure 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 
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%  varargin  command  line  arguments  to  application  (see  VARARGIN) 

O, 

o 

%  This  function  has  no  output  args,  see  OutputFcn . 

O, 

o 

%  Choose  default  command  line  output  for  application 
handles . output  =  hObject; 

%  Update  handles  structure 
guidata (  hObject,  handles  ); 

if (  strcmp (  get(hObject,  'Visible'  ),  'off'  )  ) 

initialize  gui (  hObject,  handles  ); 

end 

%  UIWAIT  makes  application  wait  for  user  response  (see  UIRESUME) 

%  uiwait (handles . figurel )  ; 

O _ 

O - 

o, 

o 

%  outputs  from  this  function  are  returned  to  the  command  line. 

O, 

o 


function  varargout  =  highResNLS  OutputFcn (  hObject,  eventdata,  handles  ) 


%  varargout 
%  hObject 
%  eventdata 
%  handles 


cell  array  for  returning  output  args  (see  VARARGOUT); 
handle  to  figure 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 


%  Get  default  command  line  output  from  handles  structure 


varargout{l}  =  handles . output; 


initialize  graphical  interface 


function  initialize  gui (  fig  handle,  handles  ) 


O, _ 

o 

o, 

o 

%  media  parameters 

o, 

o 


linearLoss  = 

0. 

0; 

gvd2deriv  = 

2  . 

2* 

(10' 

'  (-29) ) ; 

gvd3deriv  = 

0. 

0; 

linearRef ract  = 

1 . 

0; 

nonlinRef ract  = 

10 

A  ( 

-23) 

r 

%  nonlinRef ract 

= 

3* 

(10' 

'(-23));  %  from  Mechain 
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ramanSlope  =  2.0e-016; 

diffraction  =  1; 


set (  handles . linearLoss ,  'Value',  linearLoss  ); 
set (  handles . gvd2deriv,  'Value',  gvd2deriv  ); 
set (  handles . gvd3deriv,  'Value',  gvd3deriv  ); 
set (  handles . linearRefract,  'Value',  linearRef ract  ); 
set (  handles . nonlinRefract,  'Value',  nonlinRef ract  ); 
set (  handles . ramanSlope,  'Value',  ramanSlope  ); 
set (  handles . diffraction,  'Value',  diffraction  ); 

set (  handles . linearLoss ,  'String',  linearLoss  ); 
set (  handles . gvd2deriv,  'String',  gvd2deriv  ); 
set (  handles . gvd3deriv,  'String',  gvd3deriv  ); 
set (  handles . linearRefract,  'String',  linearRefract  ); 
set (  handles . nonlinRefract,  'String',  nonlinRefract  ); 
set (  handles . ramanSlope,  'String',  ramanSlope  ); 


g, _ 

o 

g, 

o 

%  beam  parameters 

g, 

o 

set (  handles . powerPaul ,  'Value',  0  ); 

set (  handles .powerSprangle,  'Value',  1  ); 
set (  handles .powerMechain,  'Value',  0  ); 

chirp  =  0.0; 

pulseEnergy  =  1 . 5*  (10A  (-2) )  ; 
wavelength  =  805 . 0* (10A (-9)  )  ; 


set (  handles . chirp,  'Value',  chirp  ); 
set (  handles . pulseEnergy,  'Value',  pulseEnergy  ); 
set (  handles . wavelength,  'Value',  wavelength  ); 

set (  handles . chirp,  'String',  chirp  ); 
set (  handles . pulseEnergy,  'String',  pulseEnergy  ); 
set (  handles . wavelength,  'String',  wavelength  ); 


g, _ 

o 

g, 

o 

%  time-axis  parameters 

g, 

o 

set (  handles . timeGrid512 ,  'Value',  1  ); 
set (  handles . timeGridl024 ,  'Value',  0  ); 

taup  =  1.25* (10A  (-14) )  ; 

widths  =  25.0; 

set (  handles . taup,  'Value',  taup  ); 
set (  handles . widths ,  'Value',  widths  ); 

set (  handles . taup,  'String',  taup  ); 
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set (  handles . widths ,  'String',  widths  ); 


O, _ 

o 

o, 

o 

%  radial-axis  parameters 

o, 

o 

set (  handles . radialGrid256,  'Value',  1  ); 

set (  handles . radialGrid512 ,  'Value',  0  ); 

set (  handles . derivFlagRawR,  'Value',  0  ); 

set (  handles . derivFlagQuarticR,  'Value',  0  ); 
set (  handles . derivFlagQuinticR,  'Value',  1  ); 

tauR  =  2.5*(10A(-3)); 

RWidth  =  25.0; 

set  (  handles . tauR,  'Value',  tauR  ); 
set (  handles . RWidth,  'Value',  RWidth  ); 

set (  handles . tauR,  'String',  tauR  ); 
set (  handles . RWidth,  'String',  RWidth  ); 

O, _ 

o 

o, 

o 

%  propagation  axis  parameters 

O, 

o 

maxPropDist  =  5.0; 
dz  =  0.00025; 

NZSlice  =  50; 


set  ( 

handles .maxPropDist, 

' Value ' , 

maxPropDist 

)  ; 

set  ( 

handles . dz , 

' Value ' , 

dz 

)  ; 

set  ( 

handles .NZSlice, 

' Value ' , 

NZSlice 

)  ; 

set  ( 

handles .maxPropDist, 

'  String' , 

maxPropDist 

)  ; 

set  ( 

handles . dz , 

'  String' , 

dz 

)  ; 

set  ( 

handles .NZSlice, 

'  String' , 

NZSlice 

) ; 

o, _ 

o 

o, 

o 

%  viewing  option 

O, 

o 

zoomOption  =  1; 

set (  handles . zoomOption,  'Value',  zoomOption  ); 

o, _ 

o 

o, 

o 

%  show  power  estimates 

o, 

o 


showPower (  handles  )  ; 
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show  power  estimates 


function  showPower (  handles  ) 


%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 


O, _ 

o 

o, 

o 

%  ref.  parameters 

O, 

o 


linearRef ract 

=  get  ( 

handles . linearRefract, 

' Value ' 

)  ; 

nonlinRef ract 

=  get  ( 

handles . nonlinRefract, 

' Value ' 

)  ; 

wavelength 

=  get  ( 

handles .wavelength. 

' Value ' 

)  ; 

pulseEnergy 

=  get  ( 

handles .pulseEnergy, 

' Value ' 

)  ; 

taup 

=  get  ( 

handles . taup. 

' Value ' 

)  ; 

tauR 

=  get  ( 

handles . tauR, 

' Value ' 

) ; 

o, _ 

o 

o, 

o 

%  peak  intensity 

o, 

o 


denom 


=  piA (3/2) * (tauRA2) *taup*erf (  1.0  )*(  1.0  -  exp (  -1.0  )  ); 


peaklntens 


=  pulseEnergy/denom; 


O, _ 

o 

o, 

o 

%  maximum  power 

O, 

o 


Pmax  =  pulseEnergy/ (2*taup) ; 

text  =  sprintf (  '  maximum  power  is  %10.3e  Watt',  Pmax  ); 

set (  handles .maxPower,  'String',  text  ); 


O, _ 

o 

o, 

o 

%  Paul  formula  for  critical  power 

O, 

o 


if (  get (  handles . powerPaul ,  'Value'  )  >  0  ) 

Pcrit  =  pi* ( ( 1 . 22 *wave length) A2 ) / ( 32 . 0*nonlinRef ract*linearRef ract ) ; 

text  =  sprintf (  '  critical  power  is  %10.3e  Watt',  Pcrit  ); 

set (  handles . criticalPower,  'String',  text  )  ; 
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%  Sprangle  formula  for  critical  power 

O, 

o 

elseif (  get  (  handles .powerSprangle,  'Value'  )  >  0  ) 

Pcrit  =  (wavelengthA2) / (2*pi*nonlinRefract*linearRefract)  ; 

text  =  sprintf (  '  critical  power  is  %10.3e  Watt',  Pcrit  ); 

set (  handles . criticalPower,  'String',  text  ); 

o, _ 

o 

o, 

o 

%  Mechain  formula  for  critical  power 

O, 

o 

else 

Pcrit  =  ( (3 . 37*wavelength) A2) / ( 8*pi*nonlinRefract*linearRefract) ; 
text  =  sprintf (  '  critical  power  is  %10.3e  Watt',  Pcrit  ); 
set (  handles . criticalPower ,  'String',  text  ); 

end 

O _ 

O - 

%  total  number  of  time  points 

O _ 

O - 

function  timeGrid512  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  timeGrid512  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'temporal  grid  size  512  selected'  ); 

set (  handles . timeGridl024 ,  'Value',  0  ); 

o _ 

O - 

function  timeGridl024^Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  timeGridl024  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'temporal  grid  size  1024  selected'  ); 

set (  handles . timeGrid512 ,  'Value',  0  ); 

o _ 

O - 

%  temporal  pulse  width 

O _ 

O - 

function  taupCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 
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%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set  (  hObject,  ' ForegroundColor '  ,  'green'  ); 

O, _ 

o 


function  taupCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

taup  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  taup  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'temporal  pulse  width  must  be  a  number',  'Error'  ); 
else 

set (  handles . taup,  'Value',  taup  ); 

end 

o _ 

O - 

%  total  number  of  pulse  widths  (time) 

O _ 

O - 


function  widthsCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  'BackgroundColor',  'black'  ); 
set (  hObject,  'ForegroundColor',  'green'  ); 

O, _ 

o 


function  widthsCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

widths  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  widths  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'number  of  pulse  widths  (time)  must  be  a  number',  'Error'  ); 
else 

set (  handles . widths ,  'Value',  widths  ); 

end 
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raw  derivative  option 


function  derivFlagRawRCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  derivFlagRawT  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'raw  (spectral)  derivative  selected  for  radial  direction'  ); 

set (  handles . derivFlagQuarticR,  'Value',  0  ); 
set (  handles . derivFlagQuinticR,  'Value',  0  ); 

o _ 

O - 

%  "quartic"  derivative  option 

O _ 

O - 


function  derivFlagQuarticR  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  derivFlagQuarticR  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

set (  handles . derivFlagRawR,  'Value',  0  ); 

set (  handles . derivFlagQuinticR,  'Value',  0  ); 

o _ 

O - 

%  "quintic"  derivative  option 

O _ 

O - 


function  derivFlagQuinticR^Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  derivFlagQuinticR  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'circularly-weighted  derivative  selected  in  radial  direction'  ); 

set (  handles . derivFlagRawR,  'Value',  0  ); 

set (  handles . derivFlagQuarticR,  'Value',  0  ); 

o _ 

O - 

%  total  number  of  radial  points 

O _ 

O - 


function  radialGrid256_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  radialGrid256  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'radial  grid  size  256  selected'  ); 


87 


set (  handles . radialGrid512 ,  'Value',  0  ); 


function  radialGrid512  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  radialGrid512  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'radial  grid  size  512  selected'  ); 


set (  handles . radialGrid256,  'Value',  0  ); 


o _ 

O - 

%  radial  pulse  width 

o _ 

O - 


function  taurCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

O, _ 

o 


function  taurCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

tauR  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  tauR  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'radial  pulse  width  must  be  a  number',  'Error'  ); 

else 

set (  handles . tauR,  'Value',  tauR  ); 
showPower (  handles  )  ; 

end 

o _ 

O - 

%  total  number  of  radial  pulse  widths 

O _ 

O - 

function  RWidthCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
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handles 


empty  -  handles  not  created  until  after  all  CreateFcns  called 


set (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 


function  RWidthCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

RWidth  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  RWidth  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'number  of  radial  pulse  widths  must  be  a  number',  'Error'  ); 
else 

set (  handles . RWidth,  'Value',  RWidth  ); 

end 

O _ 

O - 

%  maximim  propagation  distance 

O _ 

O - 


function  maxPropDistCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set (  hObject,  'BackgroundColor',  'black'  ); 
set (  hObject,  'ForegroundColor',  'green'  ); 

O, _ 

o 


function  maxPropDistCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

maxPropDist  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  maxPropDist  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'maximim  propagation  distance  must  be  a  number',  'Error'  ); 
else 

set (  handles .maxPropDist,  'Value',  maxPropDist  ); 

end 
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z-axis  step  size 


o, 

o 


function  dzCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

O, _ 

o 


function  dzCallback(  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

dz  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  dz  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'z-axis  step  size  must  be  a  number',  'Error'  ); 
else 

set (  handles. dz,  'Value',  dz  ); 

end 

o _ 

O - 

%  number  z-axis  slice  figures 

O _ 

O - 


function  NZSliceCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set (  hObject,  'BackgroundColor',  'black'  ); 
set (  hObject,  'ForegroundColor',  'green'  ); 

O, _ 

o 


function  NZSliceCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

NZSlice  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  NZSlice  )  ) 
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set (  hObject,  'String',  0  ); 


errordlg (  'number  z-axis  slice  figures  must  be  a  number',  'Error'  ); 
else 

set (  handles .NZSlice,  'Value',  NZSlice  ); 

end 

o _ 

O - 

%  linear  loss  per  meter 

O _ 

O - 


function  linearLossCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

O, _ 

o 


function  linearLossCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

linearLoss  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  linearLoss  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'linear  loss  must  be  a  number',  'Error'  ); 
else 

set (  handles . linearLoss ,  'Value',  linearLoss  ); 

end 

o _ 

O - 

%  on/off  toggler  of  2nd-order  group  vel.  deriv. 

O _ 

O - 


function  gvd2onof fCallback (  hObject,  eventdata,  handles  ) 


%  hObject 
%  eventdata 
%  handles 


handle  to  checkbox2  (see  GCBO) 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 


if (  get (  hObject,  'Value'  )  >  0  ) 

disp(  '2nd-order  group  vel.  deriv.  on 

%  gvd2deriv  =  2 . 0* (10A (-28) )  ; 

%  gvd2deriv  =  2.2*(10A(-29)); 

%  gvd2deriv  =  0 . 7*  (10A (-29) ) ; 


)  ; 

from  Schwarz 
from  Sprangle 

best  for  "leap-frog"  deriv.  (low  res.) 
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%  gvd2deriv  =  0 . 2075*  ( 10 A  ( -2 9) ) ;  %  best  for  raw  deriv.  (low  res.) 

%  gvd2deriv  =  0 . 125*  (10A  (-29) ) ;  %  too-low  value 

gvd2deriv  =  2.2*(10A(-29)); 

else 

disp(  '2nd-order  group  vel.  deriv.  off'  ); 
gvd2deriv  =  0.0; 

end 

set (  handles . gvd2deriv,  'Value',  gvd2deriv  ); 
set (  handles . gvd2deriv,  'String',  gvd2deriv  ); 

O _ 

O - 

%  2nd-order  group  vel.  deriv. 

O _ 

O - 


function  gvd2derivCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set (  hObject,  ' BackgroundColor ' ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

O, _ 

o 


function  gvd2derivCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

gvd2deriv  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  gvd2deriv  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  '2nd-order  group  vel.  deriv.  must  be  a  number',  'Error'  ); 
else 

set (  handles . gvd2deriv,  'Value',  gvd2deriv  ); 

end 

O _ 

O - 

%  3rd-order  group  vel.  deriv. 

O _ 

O - 


function  gvd3derivCreateFcn (  hObject,  eventdata,  handles  ) 


%  hObject 
%  eventdata 
%  handles 


handle  (see  GCBO) 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
empty  -  handles  not  created  until  after  all  CreateFcns  called 
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set (  hObject,  ' BackgroundColor ' , 
set (  hObject,  ' ForegroundColor ' , 


black '  ) ; 

green '  ) ; 


O, 

o 


function  gvd3derivCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

gvd3deriv  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  gvd3deriv  )  ) 

set (  hObject,  'String',  0  ); 

errordlg (  ' 3rd-order  group  vel.  deriv.  must  be  a  number',  'Error'  ); 

else 

set (  handles . gvd3deriv,  'Value',  gvd3deriv  ); 

end 

o _ 

O - 

%  linear  refractive  index 

O _ 

O - 


function  linearRef ract  CreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  linearRef ract  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

set  (  hObject,  'BackgroundColor',  'black'  ); 
set (  hObject,  'ForegroundColor',  'green'  ); 

O, _ 

o 


function  linearRef ract  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  linearRef ract  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

linearRef ract  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  linearRef ract  )  ) 

set (  hObject,  'String',  0  ); 

errordlg (  'linear  refractive  index  must  be  a  number',  'Error'  ); 
else 

set (  handles . linearRefract,  'Value',  linearRef ract  ); 
showPower (  handles  )  ; 

end 
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nonlinear  refractive  index 


function  nonlinRef ractCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

O, _ 

o 


function  nonlinRef ractCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

nonlinRef ract  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  nonlinRef ract  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'nonlinear  refractive  index  must  be  a  number',  'Error'  ); 
else 

set (  handles . nonlinRefract,  'Value',  nonlinRef ract  ); 
showPower (  handles  )  ; 

end 

o _ 

O - 

%  on/off  toggler  of  Raman  gain 

O _ 

O - 


function  ramanOnOf fCallback (  hObject,  eventdata,  handles  ) 


%  hObject 
%  eventdata 
%  handles 


handle  to  checkbox2  (see  GCBO) 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 


if (  get (  hObject,  'Value'  )  >  0  ) 
disp(  'Raman  on'  ); 


%  ramanSlope  =  7.1429e-014;  %  from  Sprangle 

%  ramanSlope  =  3.4507e-016;  %  from  Zemyanov 


ramanSlope  =  2.0e-016; 


else 

disp(  'Raman  off'  ); 


ramanSlope  =  0.0; 
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end 


set (  handles . ramanSlope,  'String',  ramanSlope  ); 
set (  handles . ramanSlope,  'Value',  ramanSlope  ); 


o, 

o 

o, 

o 

o, 

o 


slope  of  the  Raman  gain 


function  ramanSlopeCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

O, _ 

o 


function  ramanSlopeCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

ramanSlope  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  ramanSlope  )  ) 

set (  hObject,  'String',  0  ); 

errordlg (  'Raman  slope  must  be  a  number',  'Error'  ); 
else 

set (  handles . ramanSlope,  'Value',  ramanSlope  ); 

end 

o _ 

O - 

%  diffraction  on/off  selection 

O _ 

O - 


function  dif f raction_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  checkbox2  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

if (  get (  hObject,  'Value'  )  >  0  ) 
disp(  'diffraction  on'  ); 
else 

disp(  'diffraction  off'  ); 

end 

O _ 

O - 

%  chirp 

O _ 

O - 
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function  chirpCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

o, _ 

o 


function  chirpCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

chirp  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  chirp  )  ) 

set  (  hObject,  'String',  0  ); 

errordlg (  'chirp  must  be  a  number',  'Error'  ); 
else 

set (  handles . chirp,  'Value',  chirp  ); 

end 

O _ 

O - 

%  pulse  energy 

O _ 

O - 


function  pulseEnergyCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  'BackgroundColor',  'black'  ); 
set  (  hObject,  'ForegroundColor',  'green'  ); 

O, _ 

o 


function  pulseEnergyCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

pulseEnergy  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  pulseEnergy  )  ) 

set (  hObject,  'String',  0  ); 

errordlg (  'pulse  energy  must  be  a  number',  'Error'  ); 
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else 

set (  handles . pulseEnergy,  'Value',  pulseEnergy  ); 
showPower (  handles  )  ; 

end 

o _ 

O - 

%  wavelength  for  vacuum 

O _ 

O - 


function  wavelengthCreateFcn (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  empty  -  handles  not  created  until  after  all  CreateFcns  called 

set  (  hObject,  ' BackgroundColor '  ,  'black'  ); 
set (  hObject,  ' ForegroundColor '  ,  'green'  ); 

o, _ 

o 


function  wavelengthCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

wavelength  =  str2double (  get (  hObject,  'String'  )  ); 

if (  isnan (  wavelength  )  ) 

set (  hObject,  'String',  0  ); 

errordlg (  'wavelength  must  be  a  number', 'Error'  ); 
else 

set (  handles . wavelength,  'Value',  wavelength  ); 
showPower (  handles  )  ; 

end 

o _ 

O - 

%  Paul's  formulal  self-f ocusing  power 

O _ 

O - 


function  powerPaulCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  radio  button  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  ' Paul ' ' s  formula  self-f ocusing  power  option  selected'  ); 

set (  handles .powerSprangle,  'Value',  0  ); 
set (  handles .powerMechain,  'Value',  0  ); 
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showPower (  handles  )  ; 


Sprangle  formula  self-f ocusing  power 


function  powerSprangleCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  radio  button  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'Sprangle  formula  self-f ocusing  power  option  selected'  ); 

set (  handles . powerPaul ,  'Value',  0  ); 

set (  handles .powerMechain,  'Value',  0  ); 

showPower (  handles  )  ; 


Mechain  formula  (Teramobile)  self-f ocusing  power 


function  powerMechainCallback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  radio  button  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'Mechain  formula  self-f ocusing  power  option  selected'  ); 

set (  handles . powerPaul ,  'Value',  0  ); 

set (  handles .powerSprangle,  'Value',  0  ); 

gvd2deriv  =  get (  handles . gvd2deriv,  'Value'  ); 

showPower (  handles  )  ; 


o, 

o 


zoom  option 


function  zoomOption  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  checkbox2  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

if (  get (  hObject,  'Value'  )  >  0  ) 
disp (  ' zoom  on '  ) ; 

else 

disp(  'zoom  off'  ); 

end 


o, 

o 

o, 

o 


distributed  computing  option 
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function  parallel  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  parallel  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

if (  get (  hObject,  'Value'  )  >  0  ) 

disp(  'distributed  computing  on'  ); 
else 

disp(  'distributed  computing  off'  ); 

end 

O _ 

O - 

%  respond  to  button  press  (calculate) 

o _ 

O - 


function  calculateCallback (  hObject,  eventdata,  handles  ) 


%  hObject 
%  eventdata 
%  handles 


handle  to  pushbutton2  (see  GCBO) 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 


set  ( 

handles .plot3. 

' BackgroundColor 

,  'red' 

)  ; 

set  ( 

handles .plot3. 

' Enable ' , 

'off' 

)  ; 

set  ( 

handles . plot2 , 

' BackgroundColor 

,  'red' 

)  ; 

set  ( 

handles . plot2 , 

' Enable ' , 

'off' 

)  ; 

set  ( 

handles . plotl , 

' BackgroundColor 

,  'red' 

)  ; 

set  ( 

handles .plotl , 

' Enable ' , 

'off' 

)  ; 

set  ( 

handles . calculate. 

' BackgroundColor 

,  'red' 

)  ; 

set  ( 

handles . calculate. 

' Enable ' , 

'off' 

)  ; 

pause  (  0.1  )  ; 

diffraction  =  get  (  handles . diffraction. 

Value ' 

; 

if  ( 

diffraction  >  0  ) 

highResNLSfullCalc (  handles  ) ; 


else 

highResNLSf astCalc (  handles  )  ; 


end 

set  ( 

handles . calculate. 

' BackgroundColor '  , 

green ' 

) ; 

set  ( 

handles . calculate. 

' Enable ' , 

on ' 

) ; 

set  ( 

handles .plotl , 

' BackgroundColor '  , 

green ' 

) ; 

set  ( 

handles .plotl , 

' Enable ' , 

on ' 

) ; 

set  ( 

handles . plot2 , 

' BackgroundColor ' , 

green ' 

) ; 

set  ( 

handles . plot2 , 

' Enable ' , 

on ' 

) ; 

set  ( 

handles . plot3 , 

' BackgroundColor ' , 

green ' 

) ; 

set  ( 

o _ 

handles . plot3 , 

' Enable ' , 

on ' 

) ; 

o, 

o 

o _ 

respond  to  button  press 

(animation) 

99 


function  plotl  Callback (  hObject,  eventdata,  handles  ) 


%  hObject 
%  eventdata 
%  handles 


handle  to  plotl  (see  GCBO) 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 


set  ( 

handles . plot3 , 

' BackgroundColor 

r 

red '  ) ; 

set  ( 

handles .plot3. 

' Enable ' , 

off'  ); 

set  ( 

handles . plot2 , 

' BackgroundColor 

r 

red '  ) ; 

set  ( 

handles . plot2 , 

' Enable ' , 

off'  ); 

set  ( 

handles .plotl , 

' BackgroundColor 

r 

red '  ) ; 

set  ( 

handles .plotl , 

' Enable ' , 

off'  ); 

set  ( 

handles . calculate. 

' BackgroundColor 

r 

red '  ) ; 

set  ( 

handles . calculate. 

' Enable ' , 

off'  ); 

pause  (  0.1  )  ; 

highResNLSplotl ; 

set  ( 

handles . calculate. 

' BackgroundColor 

r 

green '  ) ; 

set  ( 

handles . calculate. 

' Enable ' , 

on '  )  ; 

set  ( 

handles .plotl , 

' BackgroundColor 

r 

green '  ) ; 

set  ( 

handles .plotl , 

' Enable ' , 

on '  )  ; 

set  ( 

handles . plot2 , 

' BackgroundColor 

r 

green '  ) ; 

set  ( 

handles . plot2 , 

' Enable ' , 

on '  )  ; 

set  ( 

handles . plot3 , 

' BackgroundColor 

r 

green '  ) ; 

set  ( 

handles . plot3 , 

' Enable ' , 

on '  )  ; 

o, 

o 

o _ 

respond  to  button  press 

(cumm.  plotl) 

O 


function  plot2  Callback (  hObject,  eventdata,  handles  ) 


%  hObject 
%  eventdata 
%  handles 


handle  to  plotl  (see  GCBO) 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 


set (  handles . plot3 , 
set (  handles . plot3 , 
set (  handles . plot2 , 
set (  handles . plot2 , 
set (  handles . plotl , 
set (  handles . plotl , 
set (  handles . calculate, 
set (  handles . calculate, 

pause (  0.1  ) ; 

highResNLSplot2  ; 

set  (  handles . calculate, 
set (  handles . calculate, 
set (  handles . plotl , 
set (  handles . plotl , 
set (  handles . plot2 , 


BackgroundColor'  , 

'  red ' 

)  ; 

Enable ' , 

'off' 

)  ; 

BackgroundColor' , 

'  red ' 

)  ; 

Enable ' , 

'off' 

)  ; 

BackgroundColor' , 

'  red ' 

)  ; 

Enable ' , 

'off' 

)  ; 

BackgroundColor' , 

'  red ' 

)  ; 

Enable ' , 

'off' 

) ; 

BackgroundColor' , 

' green ' 

)  ; 

Enable ' , 

'  on  ' 

)  ; 

BackgroundColor' , 

' green ' 

)  ; 

Enable ' , 

'  on  ' 

)  ; 

BackgroundColor' , 

' green ' 

) ; 
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set (  handles . plot2 , 
set (  handles . plot3 , 
set (  handles . plot3 , 


Enable ' , 

BackgroundColor ' , 
Enable ' , 


on 

green 
on  ' 


O, 

o 

o, 

o 

o, 

o 


respond  to  button  press  (cumm.  plotl) 


function  plot3_Callback (  hObject,  eventdata,  handles  ) 


%  hObject 
%  eventdata 
%  handles 


handle  to  plotl  (see  GCBO) 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 


set (  handles . plot3 , 
set (  handles . plot3 , 
set (  handles . plot2 , 
set (  handles . plot2 , 
set (  handles . plotl , 
set (  handles . plotl , 
set (  handles . calculate, 
set (  handles . calculate. 


BackgroundColor' , 

'  red 

Enable ' , 

'off 

BackgroundColor' , 

'  red 

Enable ' , 

'off 

BackgroundColor' , 

'  red 

Enable ' , 

'off 

BackgroundColor' , 

'  red 

Enable ' , 

'off 

pause (  0.1  ) ; 
highResNLSplot3 ; 


set  ( 

handles . calculate. 

' BackgroundColor ' , 

' green ' 

) ; 

set  ( 

handles . calculate. 

' Enable ' , 

'  on ' 

) ; 

set  ( 

handles .plotl , 

' BackgroundColor ' , 

' green ' 

) ; 

set  ( 

handles .plotl , 

' Enable ' , 

'  on ' 

) ; 

set  ( 

handles . plot2 , 

' BackgroundColor ' , 

' green ' 

) ; 

set  ( 

handles . plot2 , 

' Enable '  , 

'  on  ’ 

) ; 

set  ( 

handles . plot3 , 

' BackgroundColor ' , 

' green ' 

) ; 

set  ( 

o _ 

handles . plot3 , 

' Enable ' , 

'  on ' 

) ; 

o, 

o 

o _ 

end 

o 
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B.2.2  The  main  function  for  the  “fast  model”  (highResNLSfastCalc) 

For  the  “fast  model”,  the  main  program  references  highResNLSfastCalc.m: 


function  highResNLSfastCalc (  handles  ) 

S'i^^^^^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

O, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  fast  calculations  without  diffraction 

O, 

o 

%  14  Sep  2006  version  5.5 

%  with  option  for 

%  MATLAB  7  Distributed  Computing 

O, 

o 

o, 

o 

%  handles  structure  with  handles  and  user  data 

o, 

o 

%'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


format  long; 


initialize  timing  statistics 


tic; 

startCPU  =  cputime; 

O, _ 

o 


control  parameters 


%  speed  of  light  in  units  of  m/sec 

O, 

o 


speedLight  =  2 . 99792456* (10A8)  ; 


o, _ 

o 

o, 

o 

%  media  parameters 

o, 

o 


linearLoss 
gvd2deriv 
gvd3deriv 
linearRef ract 
nonlinRef ract 
ramanSlope 


get (  handles . linearLoss ,  'Value'  ); 
get (  handles . gvd2deriv,  'Value'  ); 
get (  handles . gvd3deriv,  'Value'  ); 
get (  handles . linearRefract,  'Value'  ); 
get (  handles . nonlinRefract,  'Value'  ); 
get (  handles . ramanSlope,  'Value'  ); 


disp(  'media  parameters:'  ); 
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disp  ( 

\  \  • 

sprintf (  ' 

linear  loss  term  per 

meter 

[1/m]  =  %12.6e 

)  )  r 

disp  ( 

\  ^  * 

sprintf (  ' 

2nd-order  group  vel. 

deriv . 

[secA2/m]  =  %12.6e 

)  )  r 

disp  ( 

^  ^  • 

sprintf (  ' 

3rd-order  group  vel. 

deriv . 

[secA3/m]  =  %12.6e 

)  )  r 

disp  ( 

\  \  • 

sprintf (  ' 

linear  refractive  index  (unitless)  =  %12.6e 

disp  ( 

\  \  • 

sprintf (  ' 

nonlinear  refractive 

index 

[mA2/Watt]  =  %12.6e 

disp  ( 

sprintf (  ' 

slope  of  the  Raman  gain  [sec]  =  %12.6e 

o, _ 

o 

o, 

o 

%  beam  parameters 

o, 

o 


chirp 

pulseEnergy 

wavelength 


get (  handles . chirp,  'Value'  ); 
get (  handles . pulseEnergy,  'Value'  ); 
get (  handles . wavelength,  'Value'  ); 


disp(  'beam  parameters:'  ); 

disp(  sprintf (  '  chirp  =  %4.2f'. 


disp(  '  (-1  =  down  chirp,  +1  =  up  chirp,  0 

disp(  sprintf (  '  initial  pulse  energy  [joule] 


no  chirp) '  ) ; 

=  %12 . 6e ' 


disp(  sprintf (  '  wavelength  for  vacuum  [m] 


=  %12 . 6e ' 


if (  get (  handles . powerPaul ,  'Value'  )  >  0  ) 

disp(  '  self-f ocusing  power  option  of  Paul'  ); 

elseif (  get (  handles .powerSprangle,  'Value'  )  >  0  ) 

disp(  '  self-f ocusing  power  option  of  Sprangle '  ); 

else 

disp(  '  self-f ocusing  power  option  of  Mechain '  ) ; 

end 


O, _ 

o 

o, 

o 

%  time-axis  control 

O, 

o 


if (  get (  handles . timeGridl 024 ,  'Value'  )  >  0  ) 
Ntemporal  =  1024; 
else 

Ntemporal  =  512; 

end 

taup  =  get (  handles . taup,  'Value'  ); 
widths  =  get (  handles . widths ,  'Value'  ); 

tMax  =  widths*taup; 


linearLoss 
gvd2deriv 
gvd3deriv 
linearRef ract 
nonlinRef ract 
ramanSlope 


chirp 

pulseEnergy 

wavelength 
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dt 

=  2 . 0*tMax/ (1 . 0*Ntemporal) ; 

disp  ( 

' time-axis 

parameters : '  ) ; 

disp  ( 

sprintf (  ' 

pulse  width  [sec] 

=  %12 . 6e  '  , 

taup 

) )  ; 

disp  ( 

sprintf (  ' 

total  number  of  time  points 

=  %d ' , 

Ntemporal 

) )  ; 

disp  ( 

sprintf (  ' 

total  number  of  pulse  widths 

=  %6.3f ' , 

widths 

) )  ; 

disp  ( 

sprintf (  ' 

maximum  time  [sec] 

=  %12 . 6e ' , 

tMax 

) )  ; 

disp  ( 

sprintf (  ' 

time  step  [sec] 

=  %12 . 6e  '  , 

dt 

) ) ; 

NhalfTime  = 

Ntemporal/ 2 ; 

NhalfPlus  = 

NhalfTime  +  1; 

data . 

Ntemporal  = 

Ntemporal ; 

o, _ 

o 

o, 

o 

%  radial-axis  control 

O, 

o 


if (  get (  handles . radialGrid512 ,  'Value'  )  >  0  ) 


Nradial  =  512 

r 

else 

Nradial  =  256 

r 

end 

tauR 

= 

get (  handles . tauR,  'Value'  ); 

RWidth 

= 

get (  handles . RWidth,  'Value'  ); 

Rmax 

= 

RWidth* tauR; 

NradialMinusl  = 

Nradial  -  1; 

NradialPlusl  = 

Nradial  +  1; 

dr 

= 

Rmax/ ( 1 . 0*NradialMinusl ) ; 

disp  ( 

'radial-axis  parameters:'  ); 

disp  ( 

sprintf (  ' 

FWHM  radius  [m] 

=  %12 

.  6e '  , 

tauR 

) )  ; 

disp  ( 

sprintf (  ' 

total  number  of  radial  points 

=  %d ' 

r 

Nradial 

) )  ; 

disp  ( 

sprintf (  ' 

total  number  of  radial  pulse  widths 

=  %6 . 

3f '  , 

RWidth 

) )  ; 

disp  ( 

sprintf (  ' 

maximum  radial  distance  [m] 

=  %12 

.  6e '  , 

Rmax 

) )  ; 

disp  ( 

sprintf (  ' 

radial  step  size  [m] 

=  %12 

.  6e '  , 

dr 

) ) ; 

o, _ 

o 

o, 

o 

%  propagation-axis  control 

O, 

o 


maxPropDist 

propStep 

NZSlice 


=  get (  handles .maxPropDist,  'Value'  ); 
=  get (  handles. dz,  'Value'  ); 
=  get (  handles .NZSlice,  'Value'  ); 


data . numStep 


floor (  (maxPropDist/propStep)  +  0.5  ); 


data . plotted 


f loor (  ( ( 1 . 0*data . numStep) / (1.0*NZSlice) ) 


0.5  )  ; 
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data . numPlots  =  NZSlice  +  1; 


disp(  'propagation-direction  parameters:'  ); 


disp  ( 

\  )  • 

sprintf (  ' 

maximim  propagation  distance  [m] 

=  %12 . 6e ' , 

maxPropDist 

1  )  r 

disp  ( 

\  ^  * 

sprintf (  ' 

z-axis  step  size  [m] 

=  %12 . 6e  '  , 

propStep 

)  )  r 

disp  ( 

^  ^  • 

sprintf (  ' 

total  number  of  steps  in  z 

=  %d\ 

data . numStep 

)  )  r 

disp  ( 

\  \  • 

sprintf (  ' 

number  z-axis  slice  figures 

=  %d\ 

data . numPlots 

disp  ( 

sprintf (  ' 

(plot  after  every  %dth  slice)  '  , 

data . plotted 

o, _ 

o 

o, 

o 

%  viewing  option 

O, 

o 

data . zoomOption  =  get (  handles . zoomOption,  'Value'  ); 
disp(  'plotting  parameters:'  ); 

disp(  sprintf (  '  zoom  option  ( l=zoomed, 0=wide)  =  %d',  data . zoomOption 


O, 

o 

%  temporal  plotting 

O, 

o 


if (  data . zoomOption  ==  0  ) 

data.Ntplot  =  Ntemporal/4; 
data.itplotl  =  0; 


data.itplot2  =  : 

Ntemporal ; 

else 

if (  Ntemporal  = 

=  1024  ) 

data . Ntplot 

=  257; 

data . itplotl 

=  385; 

data . itplot2 

=  641; 

else 

data . Ntplot 

=  129; 

data . itplotl 

=  193; 

data . itplot2 

=  321; 

end 

end 


o, 

o 

%  radial  plotting 

O, 

o 


if (  data . zoomOption  ==  0  ) 
Nrplot  =  Nradial/4; 
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if (  Nradial  ==  512  ) 

Nrplot  =  128; 
else 

Nrplot  =  64; 

end 

Nradial  =  Nrplot; 

NradialMinusl  =  Nradial  -  1; 

NradialPlusl  =  Nradial  +  1; 

end 

o, _ 

o 

o, 

o 

%  distribution  option 

O, 

o 

data . parallel  =  get (  handles . parallel ,  'Value'  ); 

if (  data . parallel  ==  0  ) 

disp(  '***  serial  mode  ***'  )• 

else 

disp(  '***  parallel  mode  ***'  ); 
numWorkers  =  8; 

rchunk  =  intl6(  Nradial/numWorkers  ); 

plchunk  =  intl6(  Nrplot/numWorkers  ); 

startRadius  =  1; 

endRadius  =  rchunk; 

startPlotted  =  1; 

endPlotted  =  plchunk; 

for  worker  =  1: numWorkers 

firstPlotted (worker )  =  startPlotted; 

lastPlotted (worker )  =  endPlotted; 

startPlotted  =  startPlotted  +  plchunk; 

endPlotted  =  endPlotted  +  plchunk; 

if (  data . zoomOption  >  0  ) 

firstRadius (worker)  =  firstPlotted (worker ) ; 
lastRadius (worker)  =  lastPlotted (worker ) ; 

else 

firstRadius (worker)  =  startRadius; 

lastRadius (worker)  =  endRadius; 

startRadius  =  startRadius  +  rchunk; 

endRadius  =  endRadius  +  rchunk; 

end 
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end 


end 


intensity,  power,  and  energy 


denom  =  piA (3/2) * (tauRA2) *taup*erf (  1.0  )*(  1.0  -  exp (  -1.0  )  ) 

peaklntens  =  pulseEnergy/denom; 


%  maximum  power 


Pmax 


=  pi* (tauRA2) *peaklntens; 


Paul  formula  for  critical  power 


if (  get (  handles . powerPaul ,  'Value'  )  >  0  ) 

Pcrit  =  pi* ( ( 1 . 22 *wave length) A2 ) / ( 32 . 0*nonlinRef ract*linearRef ract )  ; 


O, _ 

o 

o, 

o 

%  Sprangle  formula  for  critical  power 

O, 

o 


elseif (  get  (  handles .powerSprangle,  'Value'  )  >  0  ) 

Pcrit  =  (wavelengthA2) / (2*pi*nonlinRefract*linearRefract)  ; 


O, _ 

o 

o, 

o 

%  Mechain  formula  for  critical  power 

O, 

o 


else 

Pcrit  =  ( (3 . 37*wavelength) A2) / ( 8*pi*nonlinRefract*linearRef ract) ; 

end 


O, _ 

o 

o, 

o 

%  power  ratio  for  nonlinear  terms 

O, 

o 


powerRatio  =  Pmax/Pcrit; 
disp (  ' power : '  ) ; 

disp(  sprintf (  '  P  max  [Watt]  =  %12.6e',  Pmax  )); 
disp(  sprintf (  '  P_crit  [Watt]  =  %12.6e',  Pcrit  )); 

O, _ 

o 

o, 

o 
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%  factor  for  energy  calc. 

O, 

o 

energyFactor  =  2 . 0*pi* (peaklntens) *dr*dt; 

O, _ 

o 

%  radial  axis 

o, _ 

o 

disp(  'Forming  radial  axis  arrays'  ); 

scaledRadii  =  zeros (  Nradial,  1  ); 
trueRadii  =  zeros (  Nradial,  1  ); 

scaledDr  =  dr/tauR; 

rvalue  =  0.0; 

for  jradial  =  l:Nradial 

scaledRadii ( j radial )  =  rvalue; 

rvalue  =  rvalue  +  scaledDr; 

end 

rscaling  =  1 . 0/scaledRadii (Nradial) ; 
trueRadii  =  scaledRadii*tauR; 

o, _ 

o 

o, 

o 

%  case  A:  wide-view  plotting  axis 

o, 

o 

if (  data . zoomOption  ==  0  ) 

plotr  =  zeros  (  Nrplot,  1  ) ; 

jplot  =  1; 

for  j4  =  4:4:Nradial 

j  3  =  j4  -  1; 

j  2  =  j  3  -  1  ; 

plotr (jplot)  =  0.5* (  trueRadii ( j 2 )  +  trueRadii ( j 3 )  ); 

jplot  =  jplot  +  1; 

end 

o, _ 

o 

o, 

o 

%  case  B:  zoomed  plotting  axis 

o, 

o 

else 

plotr  =  zeros  (  Nrplot,  1  )  ; 

for  jl  =  l:Nrplot 

plotr (jl)  =  trueRadii (jl); 
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end 


end 

o, _ 

o 

%  temporal  axis 

o, _ 

o 

disp(  'Forming  temporal  axis  arrays'  ); 

scaledTimes  =  zeros (  Ntemporal,  1  ); 

trueTimes  =  zeros (  Ntemporal,  1  ); 

scaledDt  =  dt/taup; 

tvalue  =  0.0  -  (  NhalfTime  -  0.5  ) *scaledDt; 

for  itime  =  1: Ntemporal 

scaledTimes (itime)  =  tvalue; 

tvalue  =  tvalue  +  scaledDt; 

end 

trueTimes  =  scaledTimes*taup; 

o, _ 

o 

o, 

o 

%  case  A:  wide-view  plotting  axis 

o, 

o 

if (  data . zoomOption  ==  0  ) 

plott  =  zeros (  data.Ntplot,  1  ); 

i3  =  data . itplotl ; 

for  iplot  =  1 : data . Ntplot 
i3  =  i3  +  4; 

i2  =  i3  -  1; 

il  =  i2  -  1; 

plott (iplot)  =  0.5*  (  trueTimes (il)  +  trueTimes ( i2 )  ); 

end 

O, _ 

o 

o, 

o 

%  case  B:  zoomed  plotting  axis 

o, 

o 

else 

plott  =  zeros  (  data.Ntplot,  1  ); 
iplot  =  1; 

for  il  =  data . itplotl : data . itplot2 
plott (iplot)  =  trueTimes (il) ; 

iplot  =  iplot  +  1; 

end 
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end 


o, 

o 

o, 

o 


o, 

o 


pulse  energy  mask 


mask  =  zeros (  Ntemporal,  Nradial  ); 

for  jradial  =  l:Nradial 

rvalue  =  trueRadii ( j radial )  ; 

if (  rvalue  <=  tauR  ) 

for  itime  =  1: Ntemporal 

tvalue  =  abs (  trueTimes (itime)  ); 

if (  tvalue  <=  taup  ) 

mask (itime, jradial)  =  rvalue; 

end 

end 

end 

end 

o, _ 

o 

%  initial  Gaussian  pulse  definition 

O, _ 

o 

o, 

o 

%  The  pulse  is  given  in  terms  of  scaled  units  of  taup 
%  where  taup  is  the  full  width  at  half  maximum  in  intensity. 

O, 

o 

%  The  chirp  is  controlled  by  the  parameter  chirp  and  is  assumed  to  be 
%  quadratic  in  time,  i.e.,  exp (-i*chirp* (t/t  e)A2).  This  implies  that 
%  the  instantaneous  frequency  increases  linearly  from  the  leading  to 
%  trailing  edge  for  chirp  >  0  which  is  called  "up-chirp"  while  the 
%  opposite  occurs  for  chirp  <  0  which  is  called  "down-chirp" . 

O, 

o 


disp(  'Forming  Gaussian  pulse...'  ); 


O, 

o 

%  radial  part 

O, 

o 


radexpons  =  -0 . 5*scaledRadii . *scaledRadii; 
radialCurve  =  exp (  radexpons  ) ; 


o, 

o 

%  temporal  part 

O, 

o 


tratios 

timexpons 

timeCurve 


=  -0 . 5*scaledTimes . *scaledTimes; 
=  tratios* (  1.0  +  (i*chirp) ) ; 

=  exp (  timexpons  ) ; 
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full  pulse 


o, 

o 


field  =  complex ( 1 . 0 , 0 . 0 ).* (timeCurve* (radialCurve .'))  ; 


o, _ 

o 

%  derivative  operators  (time) 

O, _ 

o 


disp(  'Pre-calculating  temporal  derivative  operators...'  ); 


O, _ 

o 

o, 

o 

%  temporal  frequencies 

O, 

o 

omegas  =  zeros (  Ntemporal,  1  ); 

deltaOmega  =  2*pi/ (  dt*Ntemporal  ); 
omega  =  0.0; 

for  itime  =  l:NhalfPlus 
omegas (itime)  =  omega; 
omega  =  omega  +  deltaOmega; 

end 

itime2  =  Ntemporal; 

for  itimel  =  2:NhalfPlus 

omegas ( itime2 )  =  -omegas ( itimel ) ; 

itime2  =  itime2  -  1; 

end 

O, _ 

o 

o, 

o 

%  time-derivative  operators 

O, 

o 

data . f irstDeriv  =  i.*omegas; 

data . secondDeriv  =  omegas . *omegas ; 

data . thirdDeriv  =  i . *omegas . *omegas . *omegas ; 

o, _ 

o 

%  pre-calculated  constants 

o, _ 

o 

disp(  'Pre-calculating  constants...'  ); 

data.factorl  =  0 . 5*i*gvd2deriv*propStep; 
data.factor2  =  gvd3deriv*propStep/6 . 0 ; 

data.factor3  =  0 . 5*linearLoss*propStep; 

data.factor4  =  i* (  2*pi*nonlinRefract*propStep*peakIntens/wavelength  ); 
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data . f actor7 


data . factor4*ramanSlope; 


data.factor6  =  nonlinRefract*propStep*peakIntens/speedLight; 
data.factor5  =  2*data . factor6; 

o, _ 

o 

%  propagation  loop 

o, _ 

o 

disp(  sprintf (  'Beginning  propagaion  loop  of  %d  steps',  data.numStep  )); 

timer  =  fix(  clock  ); 

hr  =  timer ( 4 ) ; 

minv  =  timer (5); 

sec  =  timer ( 6) ; 

disp(  sprintf (  ' - >  started  at  %d:%d:%d',  hr,  minv,  sec  )); 


O, _ 

o 

o, 

o 

%  pre-allocate  arrays 

o, 

o 

intensity  =  zeros (  Nrplot,  data.Ntplot,  data . numPlots  ); 

radialEnergy  =  zeros (  Nrplot,  data . numPlots  ); 

distances  =  zeros (  data . numPlots ,  1  ); 

energies  =  zeros (  data . numPlots ,  1  ); 

o, _ 

o 

o, 

o 

%  distances 

o, 

o 

figno  =  1; 

distances (1)  =  0.0; 

for  step  =  1 : data . numStep 

propagDist  =  step*propStep; 

if (  mod (  step,  data. plotted  )  ==  0  ) 
figno  =  figno  +  1; 

distances (figno)  =  propagDist; 

end 

end 

o, _ 

o 

o, 

o 

%  set  counters 

o, 

o 


lastStep  =  data.numStep; 
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lastFrame  =  data . numPlots ; 


o, _ 

o 

o, 

o 

%  case  A:  parallel  jobs 

o, 

o 


if (  data . parallel  >  0  ) 

set (  handles . status ,  'String',  '--(distributed)--'  ); 
pause (  0.001  )  ; 


try 

for  worker 
first 
last 


1 : numWorkers 

=  firstRadius (worker) ; 
=  lastRadius (worker) ; 


data . input 
data . size 


=  f ield (:, first : last) ; 
=  last  -  first  +  1; 


data .mask 


=  mask (:, first : last); 


first 

last 


=  firstPlotted (worker ) ; 
=  lastPlotted (worker ) ; 


data . sizeplot  =  last  -  first  +  1; 
inputs {worker }  =  data; 

end 


p  j  ob 


dfevalasync ( 


QhighResNLSf astPropD, 
' LookupURL ' , 

'  FileDependencies '  , 

' StopOnError '  , 


1,  inputs, 

' WinlCM' , 

{  ' highResNLSf astPropD ' }, 
true  ) 


waitForState (  pjob,  'finished'  ); 


results  =  getAHOutputArguments  (  pjob  ); 


errmsgs  =  get (  pjob. Tasks,  { ' ErrorMessage ' }  ); 
nonempty  =  ~cellfun (  @isempty,  errmsgs  ); 
celldisp(  errmsgs (  nonempty  )  ); 

for  worker  =  1: numWorkers 

frames  =  results {worker}; 


first 

last 


=  firstPlotted (worker ) ; 
=  lastPlotted (worker); 


intensity (first : last,  :,  :) 
radialEnergy (first : last, :) 


frames . intensity; 
frames . energy; 


if (  frames . lastStep  <  lastStep  ) 
lastStep  =  frames . lastStep; 

end 
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if (  frames . lastFrame  <  lastFrame  ) 
lastFrame  =  frames . lastFrame; 

end 

end 

destroy (  p j  ob  ) ; 
catch 

data . parallel  =  0; 

disp(  '***  serial  mode  (failed  connection)  ***'  ); 

end 

end 


o, _ 

o 

o, 

o 

%  case  B:  serial  calls 

o, 

o 


if (  data . parallel  ==  0  ) 
data . input 
data . size 

data .mask 

data . sizeplot 

frames 

intensity 
radialEnergy 
lastStep 
lastFrame 

end 


o, _ 

o 

o, 

o 

%  note  singularity 

O, 

o 

if (  lastStep  <  data.numStep  ) 

disp(  ' - >  field  went  near-infinite  < - '  ); 

disp(  sprintf (  'last  step:  %d  of  %d',  lastStep,  data.numStep  )  ) 
disp(  sprintf (  'last  frame:  %d  of  %d',  lastFrame,  data . numPlots  )  ) 

end 


o, _ 

o 

o, 

o 

%  pulse  energies 

o, 

o 


=  field; 

=  Nradial; 

=  mask; 

=  Nrplot; 

=  highResNLSf astProp (  data,  handles  ); 

=  f rames . intensity; 

=  f rames . energy; 

=  frames . lastStep; 

=  frames . lastFrame; 


for  figno  =  1: lastFrame 

energies (figno)  =  trapz  (  radialEnergy (:, figno)  ) *energyFactor; 

end 

Q, 
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save  data 


o, 

o 


disp  ( 

' Saving  data 

for  output . . . '  ) ; 

save 

'  plott . mat ' 

plott; 

save 

'  plotr .mat ' 

plotr; 

save 

'  distances . mat 

'  distances; 

save 

'  energies . mat ' 

energies ; 

save 

'  intensity .mat 

'  intensity; 

o, 

o 

o, 

o 

o, _ 

display  execution  times 

disp(  'analysis  complete'  ); 
timeDiff  =  cputime  -  startCPU; 

disp(  sprintf (  '  total  client  cpu  time  =  %9.3f  SEC.',  timeDiff  )); 
disp(  sprintf (  '  total  elapsed  time  =  %9.3f  SEC.',  toe  )); 


o, 

o 

o, 

o 


end 


o, 

o 
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B.2.3  Distributed  routine  for  the  “fast  model”  (highResNLSfastPropD) 

If  the  “fast  model”  is  used  in  a  distributed  computing  environment,  the  program  uses 

highResNLSfastPropD.m : 


function  frames  =  highResNLSf astPropD (  data  ) 

S'i^^^^^'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  fast  propagation  without  diffraction  --  distributed 

O, 

o 

%  14  Sep  2006  version  5.5 

%  for 

%  MATLAB  7  Distributed  Computing 


data 


input  parameters 


%  frames  updated  stucture 

%  intensity  =  intensity  matrices  (per  figure) 

%  energy  =  pulse  energy  (per  figure) 

%  lastStep  =  last  step  calculated 

%  lastFrame  =  last  frame  number 

O, 

o 

Q-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


pre-allocate  arrays 


frames . intensity  =  zeros (  data . sizeplot,  data.Ntplot,  data . numPlots  ); 
frames . energy  =  zeros (  data . sizeplot,  data . numPlots  ); 

o, _ 


6 

o, 

o 

o, 

o 

o, 

o 

initial 

(plotted)  frame 

field 

=  data. input; 

f  igno 

=  1; 

frames 

=  fastFrame (  figno,  data. 

field,  frames  ) ; 

O, 

O - 

o, 

o 

%  propagation  loop 

o, 

o 


frames . lastStep  =  0; 

for  step  =  1 : data . numStep 
frames . lastStep  =  step; 
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for  item  =  l:data.size 

timeBuffer  =  f ield ( : , item) ; 

timeBufferSq  =  timeBuf fer . *timeBuf fer ; 

conjBuffer  =  conj (  timeBuffer  ); 

conjProd  =  con j Buffer . *timeBuf fer ; 

realpart  =  real (  timeBuffer  )  ; 

imagpart  =  imag (  timeBuffer  ); 

abslntens  =  (realpart . *realpart)  +  (imagpart .* imagpart )  ; 

intensProd  =  abslntens . *timeBuf fer; 

O, 

o 

%  forward  FFT 

O, 

o 

transField  =  fft (  timeBuffer  ) ; 

transConj  =  fft (  conjBuffer  ); 

translntens  =  fft (  abslntens  ); 

O, 

o 

%  applu  derivative  operators 

g, 

o 

firstDerivFr  =  data . firstDeriv . *transField; 
secondDerivFr  =  data . secondDeriv . *transField; 
thirdDerivFr  =  data . thirdDeriv . *transField; 

conjDerivFr  =  data . firstDeriv . *transConj ; 

intensDerivFr  =  data. firstDeriv. *translntens; 

O, 

O 

%  reverse  FFT 

g, 

o 

firstDeriv  =  ifft (  firstDerivFr  )  ; 

secondDeriv  =  ifft (  secondDerivFr  )  ; 

thirdDeriv  =  ifft (  thirdDerivFr  )  ; 

conjDeriv  =  ifft (  conjDerivFr  ); 

intensDeriv  =  ifft (  intensDerivFr  )  ; 

g, 

o 

%  change 

g, 

o 

parti  =  data . factorl . *secondDeriv; 

part2  =  data . f actor2 . *thirdDeriv; 

part3  =  data . factor3 . *timeBuf fer ; 

part4  =  data . factor4 . *intensProd; 
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part5  =  data . factor5 . *conj Prod . *firstDeriv; 

part6  =  data . f actor6 . *timeBuf ferSq . *con j Deriv; 

part7  =  data . factor7 . *timeBuf fer . *intensDeriv; 

deltaQ  =  parti  -  part2  -  part3  .  .  . 

+  part4  -  part5  -  part6  -  part7; 

f ield ( : , item)  =  timeBuffer  +  deltaQ; 

end 


o, 

o 

%  check  for  overflow 

O, 

o 


test  =  sum(  sum(  isnan (  field  ))); 

if (  test  >  0  ) 
break; 

end 


o, _ 

o 

o, 

o 

%  new  (plotted)  frame 

o, 

o 


if (  mod (  step,  data. plotted  )  ==  0  ) 
figno  =  figno  +  1; 

frames  =  fastFrame (  figno,  data,  field,  frames  ); 


o, 

o 

%  check  for  overflow 

O, 

o 


test  =  sum(  isnan (  frames . energy  )); 

if (  test  >  0  ) 
break; 

end 

end 

end 


generate  new  plot  frame 


function  frames  =  fastFrame (  figno,  data,  field,  frames  ) 


O, 

o 


o, 

o 

figno 

figure  number 

o, 

o 

data 

input  parameters 

o, 

o 

o. 

field 

latest  field  matrix 

o 

o, 

o 

frames 

updated  stucture 

intensity  =  intensity  matrices  (per  figure) 
energy  =  energy  matrices  (per  figure) 
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lastStep 

lastFrame 


last  step  calculated 
last  frame  number 


frames . lastFrame  =  figno; 


O, _ 

o 

o, 

o 

%  get  intensity 

o, 

o 


absField  =  abs  (  field  )  ; 

sqrField  =  absField . *absField; 


O, _ 

o 

o, 

o 

%  case  A:  wide  view 

o, 

o 


if (  data . zoomOption  ==  0  ) 
i4  =  data . itplotl ; 


iplot  : 

=  1 : data . Ntplot 

i4 

i4  + 

4; 

i3 

i4  - 

1; 

i2 

i3  - 

1; 

il 

i2  - 

1; 

jplot  = 

1; 

for  j4  : 

=  4:4: 

:  data . size 

j  3 

= 

\ — 1 

1 

■nT 

■r~i 

j  2 

= 

j  3  -  1; 

jl 

= 

j  2  -  1; 

magnit 

= 

sqrField (il. 

jD 

r 

magnit 

= 

magnit 

+ 

sqrField ( i2 , 

jD 

r 

magnit 

= 

magnit 

+ 

sqrField ( i3 , 

jD 

r 

magnit 

= 

magnit 

+ 

sqrField ( i4 , 

jD 

r 

magnit 

= 

magnit 

+ 

sqrField (il. 

j  2 ) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i2 , 

j  2 ) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i3 , 

j  2 ) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i4 , 

j  2 ) 

r 

magnit 

= 

magnit 

+ 

sqrField (il. 

j  3 ) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i2 , 

j  3 ) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i3 , 

j  3) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i4 , 

j  3) 

r 

magnit 

= 

magnit 

+ 

sqrField (il. 

j4) 

r 

magnit 

- 

magnit 

+ 

sqrField ( i2 , 

j4) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i3 , 

j  4 ) 

r 

magnit 

= 

magnit 

+ 

sqrField ( i4 , 

j  4 ) 

r 

frames . 

.  intensity (jplot, iplot, figno 

=  0 . 0625*magnit 
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jplot  =  jplot  +  1; 

end 


end 

suml  =  zeros  ( 

data . Ntemporal , 

l  )  ; 

sum2  =  zeros ( 

data . Ntemporal , 

l  )  ; 

sum3  =  zeros  ( 

data . Ntemporal , 

l  )  ; 

sum4  =  zeros  ( 

data . Ntemporal , 

l  ) ; 

jplot  =  1; 

for  j4  =  4 : 4 : data . size 

j  3 

\ — 1 

1 

■vT 

■|  i 

II 

j  2 

=  j3  -  1; 

jl 

\ — 1 

1 

CM 

-r~i 

II 

suml 

=  sqrField ( : , j 1 ) 

. *data . mask ( : ,  j 1 )  ; 

sum2 

=  sqrField ( : , j  2 ) 

.  *data .mask ( : ,  j 2 )  ; 

sum3 

=  sqrField ( : , j  3 ) 

.  * data .mask ( : ,  j  3 )  ; 

sum4 

=  sqrField ( : , j  4 ) 

.  *data . mask ( : ,  j  4 )  ; 

frames .energy (jplot, figno) 

=  trapz (  suml  ) 

... 

+  trapz  (  sum2  )  . . . 

+  trapz (  sum3  ) 

+  trapz (  sum4  ) ; 


jplot  =  jplot  +  1; 

end 


O, 

o 


%  case  B:  zoomed  view 

o, 

o 


else 

iplot  =  1; 

for  itime  =  data . itplotl : data . itplot2 
for  jplot  =  1 : data . sizeplot 

frames . intensity (jplot, iplot, figno)  =  sqrField (itime, jplot) 

end 

iplot  =  iplot  +  1; 

end 


suml  =  zeros (  data . Ntemporal ,  1  ); 
for  jl  =  1 : data . sizeplot 

suml  =  sqrField (:, j 1 ). *data . mask (:, j 1 )  ; 

frames .energy (j 1, figno)  =  trapz (  suml  ); 

end 

end 

O _ 
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end 
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B.2.4  Serial-mode  routine  for  the  “fast  model”  (highResNLSfastProp) 

If  the  “fast  model”  is  used  in  a  non-distributed  computing  environment,  the  program  uses 

highResNL  SfastProp.m : 


function  frames  =  highResNLSf astProp (  data,  handles  ) 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  fast  propagation  without  diffraction 

%  (client-only  version) 

O, 

o 

%  14  Sep  2006  version  5.5 

o, 

o 

input  parameters 

structure  with  handles  and  user  data 
updated  stucture 

intensity  =  intensity  matrices  (per  figure) 
energy  =  pulse  energy  (per  figure) 

lastStep  =  last  step  calculated 
lastFrame  =  last  frame  number 

O, 

o 

<i~'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  pre-allocate  arrays 

o, 

o 


%  data 
%  handles 

o, 

o 

%  frames 

O, 

o 

o, 

o 

o, 

o 

o, 

o 


frames . intensity  =  zeros (  data . sizeplot,  data.Ntplot,  data . numPlots  ); 
frames . energy  =  zeros (  data . sizeplot,  data . numPlots  ); 


6 

o, 

o 

o, 

o 

o, 

o 

initial 

(plotted)  frame 

field 

=  data. input; 

f  igno 

=  1; 

frames 

=  fastFrame (  figno,  data. 

field,  frames  ) ; 

O, 

O - 

o, 

o 

%  propagation  loop 

o, 

o 

frames . lastStep  =  0; 

for  step  =  1 : data . numStep 

if (  mod (  step,  10  )  ==  0  ) 

text  =  sprintf (  'step  %d  of  %d',  step,  data. numStep  ); 
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set  (  handles . status ,  'String',  text  ); 
pause (  0.0000001  ); 

end 

frames . lastStep  =  step; 

for  item  =  ltdata.size 

timeBuffer  =  f ield ( : , item) ; 

timeBufferSq  =  timeBuf fer . *timeBuf fer ; 

conjBuffer  =  conj (  timeBuffer  ); 

conjProd  =  con j Buffer . *timeBuf fer ; 

realpart  =  real (  timeBuffer  )  ; 

imagpart  =  imag (  timeBuffer  ); 

abslntens  =  (realpart . *realpart)  +  (imagpart .* imagpart ) ; 

intensProd  =  abslntens . *timeBuf fer; 

O, 

o 

%  forward  FFT 

g, 

o 

transField  =  fft (  timeBuffer  )  ; 

transConj  =  fft (  conjBuffer  ); 

translntens  =  fft (  abslntens  ); 

O, 

o 

%  apply  derivative  operators 

g, 

o 

firstDerivFr  =  data . f irstDeriv . *transField; 
secondDerivFr  =  data . secondDeriv . *transField; 
thirdDerivFr  =  data . thirdDeriv . *transField; 

conjDerivFr  =  data . f irstDeriv . *transConj ; 

intensDerivFr  =  data. firstDeriv. *translntens; 

g, 

o 

%  reverse  FFT 

O, 

O 

firstDeriv  =  ifft (  firstDerivFr  )  ; 

secondDeriv  =  ifft (  secondDerivFr  )  ; 

thirdDeriv  =  ifft (  thirdDerivFr  ) ; 

conjDeriv  =  ifft (  conjDerivFr  ); 

intensDeriv  =  ifft (  intensDerivFr  )  ; 
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change 


O, 

o 


parti 

part2 

part3 

part4 

part5 

part6 

part7 


=  data . factorl . *secondDeriv; 

=  data . f actor2 . *thirdDeriv; 

=  data . f actor3 . *timeBuf fer ; 

=  data . factor4 . *intensProd; 

=  data . fact or 5 . *con j  Prod . *f irstDeriv; 

=  data . fact or 6 . *timeBuf ferSq . *con j  Deriv; 
=  data . fact or 7 . *timeBuf fer . *intensDeriv; 


deltaQ 


=  parti  -  part2  -  part3 
+  part4  -  part5  -  part6  -  part7; 


f ield ( : , item)  =  timeBuffer  +  deltaQ; 

end 


o, 

o 

%  check  for  overflow 

O, 

o 


test  =  sum(  sum(  isnan (  field  ))); 

if (  test  >  0  ) 
break; 

end 


o, _ 

o 

o, 

o 

%  new  (plotted)  frame 

o, 

o 


if (  mod (  step,  data. plotted  )  ==  0  ) 
figno  =  figno  +  1; 

frames  =  fastFrame (  figno,  data,  field,  frames  ); 


o, 

o 

%  check  for  overflow 

O, 

o 


test  =  sum(  isnan (  frames . energy  )); 


if (  test  >  0  ) 
break; 

end 

end 

end 

o _ 

O - 

%  generate  new  plot  frame 

o _ 

O - 


function  frames  =  fastFrame (  figno,  data,  field,  frames  ) 
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o, 

o 


o, 

o 

figno 

figure  number 

o, 

o 

data 

input  parameters 

o, 

o 

o. 

field 

latest  field  matrix 

o 

o, 

o 

frames 

updated  stucture 

o, 

o 

intensity  = 

intensity  matrices  (per  figure 

o, 

o 

energy  = 

energy  matrices  (per  figure) 

o, 

o 

lastStep  = 

last  step  calculated 

o, 

o 

lastFrame  = 

last  frame  number 

o, 

o 


frames . lastFrame  =  figno; 


O, _ 

o 

o, 

o 

%  get  intensity 

o, 

o 

absField  =  abs  (  field  )  ; 

sqrField  =  absField . *absField; 

o, _ 

o 

o, 

o 

%  case  A:  wide  view 

o, 

o 


if (  data . zoomOption  ==  0  ) 


i4  =  data. 

itplotl ; 

for  iplot 

=  1 : data . Ntplot 

i4 

i4  +  4; 

i3 

i4  -  1; 

i2 

i3  -  1; 

il 

i2  -  1; 

jplot  = 

1; 

for  j  4 

=  4:4: data . size 

j  3 

=  j  4  -  1; 

j  2 

=  j  3  -  1; 

jl 

\ — 1 

1 

CM 

■r~i 

II 

magnit  =  sqrField (il, j 1) ; 
magnit  =  magnit  +  sqrField ( i2 , j 1 ) ; 
magnit  =  magnit  +  sqrField ( i3 , j 1 ) ; 
magnit  =  magnit  +  sqrField ( i4 , j 1 ) ; 

magnit  =  magnit  +  sqrField (il, j  2 ) ; 
magnit  =  magnit  +  sqrField ( i2 , j 2 ) ; 
magnit  =  magnit  +  sqrField ( i3 , j 2 ) ; 
magnit  =  magnit  +  sqrField ( i4 , j 2 ) ; 

magnit  =  magnit  +  sqrField ( il , j 3 ) ; 
magnit  =  magnit  +  sqrField ( i2 , j 3 ) ; 
magnit  =  magnit  +  sqrField ( i3 , j 3); 
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magnit  =  magnit  +  sqrField ( i4 , j 3 ) ; 

magnit  =  magnit  +  sqrField ( il , j 4 ) ; 
magnit  =  magnit  +  sqrField ( i2  ,  j 4  )  ; 
magnit  =  magnit  +  sqrField ( i3 , j 4 ) ; 
magnit  =  magnit  +  sqrField ( i4 , j 4 ) ; 

frames . intensity (jplot, iplot, figno)  =  0 . 0625*magnit; 


jplot 

end 


end 

suml  = 

zeros ( 

data . Ntemporal 

sum2  = 

zeros ( 

data . Ntemporal 

sum3  = 

zeros ( 

data . Ntemporal 

sum4  = 

zeros ( 

data . Ntemporal 

jplot  = 

1; 

for  j4  : 

=  4:4: data . size 

j  3 
j  2 
jl 


suml  = 

sum2  = 

sum3  = 

sum4  = 

frames .energy (jplot,  figno)  = 

+ 

+ 

+ 

jplot  = 


=  jplot  +  1; 


)  ; 
)  ; 
)  ; 
)  ; 


j  4  -  1; 

j  3  -  1  ; 

j  2  -  1  ; 

sqrField ( : , j 1 ) 

.  *data . mask ( : , j 1 ) ; 

sqrField ( : , j  2 ) 

.  *data .mask ( : ,  j 2 ) ; 

sqrField ( : , j  3 ) 

.  *data . mask ( : , j  3 ) ; 

sqrField ( : , j  4 ) 

.  *data .mask ( : , j  4 )  ; 

trapz (  suml  ) 
trapz  (  sum2  ) 
trapz (  sum3  ) 
trapz (  sum4  ) ; 

jplot  +  1; 

end 


o, _ 

o 

o, 

o 

%  case  B:  zoomed  view 

o, 

o 


else 

iplot  =  1; 

for  itime  =  data . itplotl : data . itplot2 
for  jplot  =  1 : data . sizeplot 

frames . intensity (jplot, iplot, figno)  =  sqrField (itime, jplot) 

end 

iplot  =  iplot  +  1; 

end 

suml  =  zeros (  data . Ntemporal ,  1  ); 
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for  jl  =  1 : data . sizeplot 
suml 


=  sqrField ( : , j 1 ) . *data . mask ( : , j 1 ) ; 


frames .energy (j 1, figno)  =  trapz (  suml  ); 
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B.2.5  The  main  function  for  the  “slower  model”  (highResNLSfullCalc) 

For  the  “slower  model”,  the  main  program  references  highResNLSfullCalc.m: 


function  highResNLSfullCalc (  handles  ) 

9'******************'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  full  calculations  with  diffraction 

O, 

o 

%  14  Sep  2006  version  5.5 

%  with  option  for 

%  MATLAB  7  Distributed  Computing 

O, 

o 

o, 

o 

%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

O, 

o 

<H-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 


format  long; 


initialize  timing  statistics 


tic; 


startCPU  =  cputime; 


O, _ 

o 

o, 

o 

o, _ 

o 

o, 

o 

%  speed  of  light  in  units 

O, 

o 


control  parameters 


of  m/sec 


speedLight  =  2 . 99792456* (10A8) ; 


o, _ 

o 

o, 

o 

%  media  parameters 

o, 

o 


linearLoss 
gvd2deriv 
gvd3deriv 
linearRef ract 
nonlinRef ract 
ramanSlope 


get (  handles . linearLoss , 
get (  handles . gvd2deriv, 
get (  handles . gvd3deriv, 
get (  handles . linearRefract, 
get (  handles . nonlinRefract, 
get (  handles . ramanSlope, 


'  Value '  ) ; 
'  Value '  ) ; 
'  Value '  ) ; 
'  Value '  ) ; 
'  Value '  ) ; 
'  Value '  ) ; 


disp(  'media  parameters:'  ); 

disp(  sprintf (  '  linear  loss  term  per  meter  [1/m]  =  %12.6e',  linearLoss 
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2nd-order  group  vel.  deriv.  [secA2/m]  =  %12.6e',  gvd2deriv 
3rd-order  group  vel.  deriv.  [secA3/m]  =  %12.6e',  gvd3deriv 
linear  refractive  index  (unitless)  =  %12.6e',  linearRef ract 
nonlinear  refractive  index  [mA2/Watt]  =  %12.6e',  nonlinRef ract 
slope  of  the  Raman  gain  [sec]  =  %12.6e',  ramanSlope 


O, _ 

o 

o, 

o 

%  beam  parameters 

o, 

o 

chirp  =  get (  handles . chirp,  'Value'  ); 

pulseEnergy  =  get (  handles . pulseEnergy,  'Value'  ); 

wavelength  =  get (  handles . wavelength,  'Value'  ); 

disp(  'beam  parameters:'  ); 
disp(  sprintf (  '  chirp 
) )  ; 

disp(  '  (-1  =  down  chirp,  +1  =  up  chirp,  0 

disp(  sprintf (  '  initial  pulse  energy  [joule] 

) )  ; 

disp(  sprintf (  '  wavelength  for  vacuum  [m] 


if (  get (  handles . powerPaul ,  'Value'  )  >  0  ) 

disp(  '  self-f ocusing  power  option  of  Paul'  ); 

elseif (  get (  handles .powerSprangle,  'Value'  )  >  0  ) 

disp(  '  self-f ocusing  power  option  of  Sprangle '  ); 

else 

disp(  '  self-f ocusing  power  option  of  Mechain '  ); 

end 

O, _ 

o 

o, 

o 

%  time-axis  control 

O, 

o 

if (  get (  handles . timeGridl 024 ,  'Value'  )  >  0  ) 
Ntemporal  =  1024; 
else 

Ntemporal  =  512; 

end 

taup  =  get (  handles . taup,  'Value'  ); 

widths  =  get (  handles . widths ,  'Value'  ); 

tMax  =  widths*taup; 

dt  =  2 . 0*tMax/ ( 1 . 0*Ntemporal ) ; 


=  %4 . 2f ' ,  chirp 

=  no  chirp) '  ) ; 

=  %12.6e',  pulseEnergy 

=  %12.6e',  wavelength 


disp (  sprintf (  ' 
) )  ; 

disp (  sprintf (  ' 
)  )  ; 

disp (  sprintf (  ' 
) )  ; 

disp (  sprintf (  ' 
) )  ; 

disp (  sprintf (  ' 
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disp(  'time-axis 

parameters : '  ) ; 

disp (  sprintf  (  ' 

pulse  width  [sec] 

=  %12 . 6e ' , 

taup 

) )  ; 

disp (  sprintf (  ' 

total  number  of  time  points 

=  %d ' , 

Ntemporal 

) )  ; 

disp (  sprintf (  ' 

total  number  of  pulse  widths 

=  %6.3f ' , 

widths 

) )  ; 

disp (  sprintf (  ' 

maximum  time  [sec] 

=  %12 . 6e ' , 

tMax 

) )  ; 

disp (  sprintf  (  ' 

time  step  [sec] 

=  %12 . 6e  '  , 

dt 

) ) ; 

Nhalf Time 

=  Ntemporal /2; 

Nhalf Plus 

=  Nhalf Time  +  1; 

dataB . Ntemporal 

=  Ntemporal; 

o, _ 

o 

o, 

o 

%  radial-axis  control 

O, 

o 


if (  get (  handles . radialGrid512 ,  'Value'  )  >  0  ) 

Nradial  =  512; 
else 

Nradial  =  256; 

end 

if (  get (  handles . derivFlagRawR,  'Value'  )  >  0  ) 

derivFlagR  =  1; 

elseif (  get  (  handles . derivFlagQuarticR,  'Value'  )  >  0  ) 
derivFlagR  =2; 
else 

derivFlagR  =  3; 

end 

tauR  =  get (  handles . tauR,  'Value'  ); 

RWidth  =  get (  handles . RWidth,  'Value'  ); 

Rmax  =  RWidth*tauR; 

NradialPlusl  =  Nradial  +  1; 


%12 . 6e ' ,  tauR  ) ) ; 
%d ' ,  Nradial  ) ) ; 
%6 . 3f ' ,  RWidth  ) ) ; 
%12 . 6e ' ,  Rmax  ) ) ; 
%d'  ,  derivFlagR  ) ) ; 

dataA. Nradial  =  Nradial; 


disp(  'radial-axis  parameters:'  ) 
disp (  sprintf ( 
disp (  sprintf ( 
disp (  sprintf ( 
disp (  sprintf ( 
disp (  sprintf ( 


FWHM  radius  [m] 
total  number  of  radial  points 
total  number  of  radial  pulse  widths 
maximum  radial  distance  [m] 
derivative  option 


O, _ 

o 

o, 

o 

%  propagation-axis  control 

O, 

o 


maxPropDist  =  get (  handles .maxPropDist,  'Value'  ); 

propStep  =  get (  handles. dz,  'Value'  ); 
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NZSlice 

=  get (  handles .NZSlice,  'Value'  ); 

numStep 

=  floor (  (maxPropDist/propStep)  +  0.5 

)  ; 

plotted 

=  floor(  ( ( 1 . 0*numStep) / ( 1 . 0*NZSlice) ) 

+  0.5  )  ; 

plotting . numPlots 

=  NZSlice  +  1; 

disp(  'propagation-direction  parameters:'  ); 
disp(  sprintf (  '  maximim  propagation  distance  [m] 

i  i  . 

=  %12 . 6e  '  , 

/  )  r 

disp (  sprintf (  ' 

\  ^  * 

z-axis  step  size  [m] 

=  %12 . 6e ' , 

)  )  r 

disp (  sprintf  (  ' 

^  ^  • 

total  number  of  steps  in  z 

=  %d'  , 

)  )  r 

disp (  sprintf (  ' 

number  z-axis  slice  figures 

=  %d'  , 

plotting . numPlots  )); 

disp(  sprintf (  '  (plot  after  every  %dth  slice) 


maxPropDist 

propStep 

numStep 

plotted 


o, _ 

o 

o, 

o 

%  viewing  option 

O, 

o 


plotting . zoomOption  =  get (  handles . zoomOption,  'Value'  ); 
disp(  'plotting  parameters:'  ); 

disp(  sprintf (  '  zoom  option  ( l=zoomed, 0=wide)  =  %d', 

plotting . zoomOption  )); 


O, 

o 

%  temporal  plotting 

O, 

o 


if (  plotting . zoomOption  ==  0  ) 

plotting . Ntplot  =  Ntemporal/4; 
plotting . itplotl  =  0; 

plotting . itplot2  =  Ntemporal; 

else 


if (  Ntemporal  ==  1024  ) 


plotting .Ntplot 

=  257; 

plotting .itplotl 
plotting . itplot2 

=  385; 
=  641; 

else 

plotting .Ntplot 

=  129; 

plotting .itplotl 
plotting . itplot2 

end 

=  193; 
=  321; 

end 


O, 

o 
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radial  plotting 


O, 

o 


if (  plotting . zoomOption  ==  0  ) 
plotting . Nrplot  =  Nradial/4; 

else 

if (  Nradial  ==  512  ) 

plotting . Nrplot  =  128; 
else 

plotting . Nrplot  =  64; 

end 

end 

O, _ 

o 

o, 

o 

%  distribution  option 

O, 

o 


parallel  =  get (  handles . parallel ,  'Value'  ); 

if (  parallel  ==  0  ) 

disp(  '***  serial  mode  ***'  ); 


else 

disp(  '***  parallel  mode  ***'  ); 


numWorkers  =  8; 


tchunk  =  intl6(  Ntemporal/numWorkers  ); 

rchunk  =  intl6(  Nradial/numWorkers  ); 


startTime  =  1; 

endTime  =  tchunk; 


startRadius  =  1; 
endRadius  =  rchunk; 


for  worker  =  1: numWorkers 

firstTime (worker)  =  startTime; 

lastTime (worker)  =  endTime; 


startTime 

endTime 

firstRadius (worker) 
lastRadius (worker) 

startRadius 

endRadius 

end 

end 


startTime  +  tchunk; 
endTime  +  tchunk; 

startRadius ; 
endRadius ; 

startRadius  +  rchunk; 
endRadius  +  rchunk; 


O, 

o 

o, 

o 

o. 


intensity,  power,  and  energy 
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denom  =  piA (3/2) * (tauRA2) *taup*erf (  1.0  )*(  1.0  -  exp (  -1.0  )  ); 

peaklntens  =  pulseEnergy/denom; 


g, 

o 


%  maximum  power 

O, 

O 


Pmax  =  pi* (tauRA2) *peaklntens; 


g. _ 

o 

o, 

o 

%  Paul  formula  for  critical  power 

Q, 

O 


if (  get (  handles . powerPaul ,  'Value'  )  >  0  ) 

Pcrit  =  pi* ( ( 1 . 22 *wave length) A2 ) / ( 32 . 0*nonlinRef ract*linearRef ract ) 


O, _ 

o 

o, 

o 

%  Sprangle  formula  for  critical  power 

g, 

o 


elseif (  get (  handles .powerSprangle,  'Value'  )  >  0  ) 

Pcrit  =  (wavelengthA2 ) / (2*pi*nonlinRefract*linearRefract)  ; 


g, _ 

o 

g, 

o 

%  Mechain  formula  for  critical  power 

g, 

o 


else 

Pcrit  =  ( (3 . 37*wavelength) A2) / ( 8*pi*nonlinRefract*linearRefract) ; 

end 


g, _ 

o 

g, 

o 

%  power  ratio  for  nonlinear  terms 

g, 

o 


powerRatio  =  Pmax/Pcrit; 


disp (  ' power : ' 
disp (  sprintf ( 
disp (  sprintf ( 


P  max  [Watt] 
P  crit  [Watt] 


%12 . 6e ' ,  Pmax  ) ) ; 
%12 . 6e ' ,  Pcrit  ) ) ; 


g, _ 

o 

g, 

o 

%  factor  for  energy  calc. 

g, 

o 


plotting . factor  =  2 . 0*pi* (peaklntens) *dt; 


o_ 
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DHT  initialization 


o 


disp(  'Pre-calculating  arrays  for  DHT...'  ); 


O, _ 

o 

o, 

o 

%  read  in  coefficients  matrix 

O, 

o 


disp(  ' 
load  c.mat; 
coef  s 
cscaling 


(loading  coefficients  matrix) ' 

=  c ( 1 , 1 : Nradial )  ; 

=  1 . 0/c  (1 , NradialPlusl )  ; 


O, _ 

o 

o, 

o 

%  radii  vector  (non-uniform  spacing) 

O, 

o 


trueRadii  =  (coef s ' ) *Rmax*cscaling; 


O, _ 

o 

o, 

o 

%  frequency  vector  (non-uniform  spacing) 

O, 

o 


rscaling  =  1.0/Rmax; 

radOmegas  =  (coef s ' ) *rscaling; 


O, _ 

o 

o, 

o 

%  XM  is  the  transformation  matrix 

O, 

o 


disp  ( 


(transformation  matrix) '  ) ; 


[ Jn, Jm] 

xdenom 

xarg 

dataA.XM 


=  meshgrid(  coefs,  coefs  ); 

=  abs (  besselj (  1,  Jn  )  ).*abs(  besselj (  1, 

=  Jn . * Jm. *cscaling; 

=  (2 . 0*cscaling) *bessel j (  0,  xarg  ). /xdenom; 


O, _ 

o 

o, 

o 

%  DHT  scalings 

O, 

o 


disp  ( 


(scaling  matrices) '  ) ; 


dataA.ml  =  (  abs (  besselj (  1,  coefs 

dataA.mlrecip  =  1 . 0 ./dataA.ml; 


)  . *rscaling  )  '  ; 


Jm  )  )  ; 
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radial  axis 


o, _ 

o 

disp(  'Forming  radial  axis  arrays'  ); 
plotting . radii  =  trueRadii; 
scaledRadii  =  trueRadii . /tauR; 

o, _ 

o 

o, 

o 

%  case  A:  wide-view  plotting  axis 

o, 

o 

if (  plotting . zoomOption  ==  0  ) 

plotr  =  zeros  (  plotting. Nrplot,  1  ); 

jplot  =  1; 

for  j4  =  4:4:Nradial 

j  3  =  j4  -  1; 

j  2  =  j  3  -  1  ; 

plotr (jplot)  =  0.5* (  trueRadii ( j 2 )  +  trueRadii ( j 3 )  ); 

jplot  =  jplot  +  1; 

end 

o, _ 

o 

o, 

o 

%  case  B:  zoomed  plotting  axis 

o, 

o 

else 

plotr  =  zeros  (  plotting. Nrplot,  1  ); 

for  jl  =  1 : plotting . Nrplot 
plotr (jl)  =  trueRadii ( j 1 ) ; 

end 

end 

O, _ 

o 

%  temporal  axis 

O, _ 

o 

disp(  'Forming  temporal  axis  arrays'  ); 

scaledTimes  =  zeros (  Ntemporal,  1  ); 

trueTimes  =  zeros (  Ntemporal,  1  ); 

scaledDt  =  dt/taup; 

tvalue  =  0.0  -  (  NhalfTime  -  0.5  ) *scaledDt; 

for  itime  =  l:Ntemporal 
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scaledTimes (itime)  =  tvalue; 

tvalue  =  tvalue  +  scaledDt; 

end 

trueTimes  =  scaledTimes*taup; 


o, _ 

o 

o, 

o 

%  case  A:  wide-view  plotting  axis 

o, 

o 


if (  plotting . zoomOption  ==  0  ) 

plott  =  zeros (  plotting. Ntplot,  1  ); 

i3  =  plotting . itplotl ; 

for  iplot  =  1 : plotting . Ntplot 
i 3  =  i  3  +  4 ; 

i2  =  i3  -  1; 

il  =  i2  -  1; 

plott (iplot)  =  0.5* (  trueTimes (il)  +  trueTimes ( i2 )  ); 

end 


O, _ 

o 

o, 

o 

%  case  B:  zoomed  plotting  axis 

o, 

o 


else 

plott  =  zeros  (  plotting. Ntplot,  1  ); 


iplot  =  1; 


for  il  =  plotting . itplotl : plotting . itplot2 
plott (iplot)  =  trueTimes (il) ; 

iplot  =  iplot  +  1; 

end 

end 

O, _ 

o 

%  pulse  energy  mask 

O, _ 

o 


plotting . Nradial  =  Nradial; 

plotting .mask  =  zeros (  Ntemporal,  Nradial  ); 

for  jradial  =  l:Nradial 

rvalue  =  trueRadii ( j radial ) ; 

if (  rvalue  <=  tauR  ) 

for  itime  =  l:Ntemporal 

tvalue  =  abs (  trueTimes (itime)  ); 
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rvalue; 


if (  tvalue  <=  taup  ) 

plotting . mask ( itime, j radial )  = 

end 

end 

end 

end 

o, _ 

o 

%  initial  Gaussian  pulse  definition 

O, _ 

o 

o, 

o 

%  The  pulse  is  given  in  terms  of  scaled  units  of  taup 
%  where  taup  is  the  full  width  at  half  maximum  in  intensity. 

O, 

o 

%  The  chirp  is  controlled  by  the  parameter  chirp  and  is  assumed  to  be 
%  quadratic  in  time,  i.e.,  exp (-i*chirp* (t/t  e)A2).  This  implies  that 
%  the  instantaneous  frequency  increases  linearly  from  the  leading  to 
%  trailing  edge  for  chirp  >  0  which  is  called  "up-chirp"  while  the 
%  opposite  occurs  for  chirp  <  0  which  is  called  "down-chirp". 

O, 

o 


disp(  'Forming  Gaussian  pulse...'  ); 


O, 

o 

%  radial  part 

O, 

o 


radexpons  =  -0 . 5*scaledRadii . *scaledRadii; 
radialCurve  =  exp (  radexpons  ) ; 


o, 

o 

%  temporal  part 

O, 

o 


tratios 

timexpons 

timeCurve 


=  -0 . 5*scaledTimes . *scaledTimes; 
=  tratios* (  1.0  +  (i*chirp) ) ; 

=  exp (  timexpons  ) ; 


O, 

o 

%  full  pulse 

o, 

o 


field  =  complex ( 1 . 0 , 0 . 0 ).* (timeCurve* ( radialCurve .')) ; 

o,  _ 

o 

%  derivative  operators  (time) 

O,  _ 

O 

disp(  'Pre-calculating  temporal  derivative  operators...'  ); 


temporal  frequencies 


137 


omegas 


zeros (  Ntemporal,  1  ); 


deltaOmega  =  2*pi/ (  dt*Ntemporal  ) ; 
omega  =  0.0; 

for  itime  =  l:NhalfPlus 
omegas (itime)  =  omega; 
omega  =  omega  +  deltaOmega; 

end 

itime2  =  Ntemporal; 

for  itimel  =  2:NhalfPlus 

omegas ( itime2 )  =  -omegas ( itimel ) ; 

itime2  =  itime2  -  1; 

end 

o, _ 

o 

o, 

o 

%  time-derivative  operators 

O, 

o 

dataB . f irstDeriv  =  i.*omegas; 

dataB . secondDeriv  =  omegas . *omegas ; 

dataB . thirdDeriv  =  i . *omegas . *omegas . *omegas ; 

o, _ 

o 

%  pre-calculated  constants 

o, _ 

o 

disp(  'Pre-calculating  temporal  constants...'  ); 

dataB . factorl  =  0 . 5*i*gvd2deriv*propStep; 
dataB . factor2  =  gvd3deriv*propStep/6 . 0 ; 

dataB . factor3  =  0 . 5*linearLoss*propStep; 

dataB . factor4  =  i* (  2*pi*nonlinRefract*propStep*peakIntens/wavelength  ); 
dataB . factor7  =  dataB . factor4*ramanSlope; 

dataB . factor6  =  nonlinRefract*propStep*peakIntens/speedLight; 
dataB . factor5  =  2*dataB . factor6; 

o, _ 

o 

%  diffraction  operator 

o, _ 

o 

disp(  'Pre-calculating  factors  for  diffraction...'  ); 


O, 

o 

o, 

o 
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radial  derivative  operator 


g, 

o 


dataA. derivOp  =  i . *radOmegas ; 


g, 

o 


%  weighting  functions  for  radial  derivative 

g, 

o 


if (  derivFlagR  ==  2  ) 

rNqu  =  radOmegas (Nradial) A4; 

denom  =  1.0 /rNqu; 

for  j radial  =  lrNradial 

rqu  =  radOmegas ( j radial ) A4 ; 

ratio  =  (  rNqu  -  rqu  ) *denom; 

dataA. derivOp ( j radial )  =  dataA. derivOp ( j radial ). *ratio; 

end 

elseif  (  derivFlagR  ==  3  ) 

rNqu  =  abs (  radOmegas (Nradial) A5  ); 

denom  =  1.0 /rNqu; 

for  j radial  =  lrNradial 

rqu  =  abs (  radOmegas (j radial) A5  ); 

ratio  =  (  rNqu  -  rqu  ) *denom; 

dataA. derivOp ( j radial )  =  dataA. derivOp ( j radial ). *ratio; 

end 

end 

O. _ 

o 

o, 

o 

%  scaling  factors  for  diffraction 

Q, 

O 


dataA . divRad 


1.0./ trueRadii  ; 


f actorV 
f actorD 


=  ( i*wavelength*propStep) / (  4*pi*linearRefract  ); 
=  factorV* (  wavelength/ (  2*pi*speedLight  ) ) ; 


dataB . scalel 
dataB . scale2 
dataB . scale3 
dataB . scale4 


factorV; 
factorV; 
f actorD; 
f actorD; 


o, 

o 

g, 

o 

o. 


initial  frame 
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disp(  'Generating  initial  (plotted)  frame...'  ); 

frames . intensity  =  zeros (  plotting. Nrplot,  ... 

plotting. Ntplot,  ... 
plotting . numPlots  ); 

frames . energy  =  zeros (  plotting . numPlots ,  1  ); 
distances  =  zeros (  plotting . numPlots ,  1  ); 

figno  =  1; 

distances (1)  =  0.0; 

frames  =  generateFrame (  figno,  plotting,  field,  frames  ); 

lastStep  =  0; 

o, _ 

o 

%  pre-allocate  arrays 

o, _ 

o 

if (  parallel  >  0  ) 

disp(  'pre-allocating  arrarys  for  distributed  processing...'  ); 

firstDeriv  =  complex ( 1 . 0 , 0 . 0 ).* zeros (  Ntemporal,  Nradial  ); 
secondDeriv  =  complex ( 1 . 0 , 0 . 0 ).* zeros (  Ntemporal,  Nradial  ); 

end 

o, _ 

o 

%  propagation  loop 

o, _ 

o 

disp(  sprintf (  'Beginning  propagaion  loop  of  %d  steps',  numStep  )); 

timer  =  fix(  clock  ); 

hr  =  timer ( 4 ) ; 

minv  =  timer (5); 

sec  =  timer ( 6)  ; 

disp(  sprintf (  ' - >  started  at  %d:%d:%d',  hr,  minv,  sec  )); 

broken  =  0; 

for  step  =  1: numStep 

text  =  sprintf (  'step  %d  of  %d',  step,  numStep  ); 

set (  handles . status ,  'String',  text  ); 
pause (  0 . 00001  ) ; 
lastStep  =  step; 

O, _ 

O 

O, 

O 

%  case  A:  parallel  jobs 
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o, 

o 

if (  parallel  >  0  ) 
try 

for  worker  =  1 : numWorkers 

first  =  firstTime (worker) ; 

last  =  lastTime (worker) ; 

dataA.size  =  last  -  first  +  1; 

dataA. input  =  field (first :  last,  :)  ; 

inputsA{ worker }  =  dataA; 

end 

pjobA  =  dfevalasync (  QhighResNLSfullPartA,  1,  inputsA,  ... 

'LookupURL',  ' WinlCM' , 

' FileDependencies ' ,  { ' highResNLSfullPartA ' } ,  ... 

' StopOnError ' ,  true  ); 

waitForState (  pjobA,  'finished'  ); 

resultsA  =  getAHOutputArguments  (  pjobA  ); 

errmsgsA  =  get (  pjobA. Tasks,  { ' ErrorMessage ' }  ); 
nonemptyA  =  ~cellfun (  @isempty,  errmsgsA  ); 
celldisp(  errmsgsA (  nonemptyA  )  ); 

for  worker  =  1: numWorkers 

returnedA  =  resultsAj worker } ; 

first  =  firstTime (worker) ; 

last  =  lastTime (worker) ; 

firstDeriv (first : last, ; )  =  returnedA. firstDeriv; 

secondDeriv (first : last, : )  =  returnedA. secondDeriv; 

end 

destroy (  pjobA  ); 

o, _ 

o 

for  worker  =  1: numWorkers 

first  =  firstRadius (worker) ; 

last  =  lastRadius (worker) ; 

dataB.size  =  last  -  first  +  1; 

dataB. input  =  f ield (:,  first : last)  ; 

dataB . f irstRadDeriv  =  f irstDeriv (:, first : last) ; 
dataB . secondRadDeriv  =  secondDeriv (:, first : last) ; 

inputsB{worker }  =  dataB; 

end 

pjobB  =  dfevalasync (  QhighResNLSfullPartB,  1,  inputsB,  ... 
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'LookupURL',  ' WinlCM' , 

' FileDependencies '  ,  { ' highResNLSfullPartB ' } , 

' StopOnError ' ,  true  ) 

waitForState (  pjobB,  'finished'  ); 

resultsB  =  getAHOutputArguments  (  pjobB  ); 

errmsgsB  =  get (  pjobB. Tasks,  { ' ErrorMessage ' }  ); 
nonemptyB  =  ~cellfun (  @isempty,  errmsgsB  ); 
celldisp(  errmsgsB (  nonemptyB  )  ); 

for  worker  =  lrnumWorkers 

first  =  firstRadius (worker) ; 

last  =  lastRadius (worker) ; 

f ield (:, first : last)  =  resultsB{worker } ; 

end 

destroy (  pjobB  ); 

o, _ 

o 


catch 

parallel  =  0; 

disp(  '***  serial  mode  (failed  connection)  ***'  )• 

end 

end 


o. _ 

o 

o, 

o 

%  case  B:  serial  calls 

o, 

o 


if (  parallel  ==  0  ) 

dataA. input  =  field; 

dataA.size  =  Ntemporal; 


returnedA  =  highResNLSfullPartA (  dataA  ); 

dataB. input  =  field; 

dataB . f irstRadDeriv  =  returnedA. firstDeriv; 
dataB . secondRadDeriv  =  returnedA. secondDeriv; 


dataB . size 


=  Nradial; 


field 

end 


=  highResNLSfullPartB (  dataB  ); 


test  =  sum(  sum(  isnan (  field  ))); 

if (  test  >  0  ) 
broken  =  1; 
break; 
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end 


o 

o, 

o 

%  new  frame 

O, 

o 

if (  mod (  step,  plotted  )  ==  0  ) 


f  igno 

=  figno  +  1; 

distances (figno) 

=  step*propStep; 

frames 

=  generateFrame (  figno,  plotting,  field,  frames  ); 

test 

=  sum(  isnan (  frames . energy  )); 

if (  test  >  0  ) 

broken  =  1; 

break; 

end 


end 

end 

o, _ 

o 

o, 

o 

%  note  singularity 

O, 

o 


if (  broken  ==  1  ) 

disp(  ' - >  field  went  near-infinite  < - '  ); 

disp(  sprintf (  'last  step:  %d  of  %d',  lastStep,  numStep  )  ); 

disp(  sprintf (  'last  frame:  %d  of  %d',  frames . lastFrame,  plotting . numPlots  ) 


end 


O, 

o 

o, 

o 


o, 

o 


save  data 


disp(  'Saving  data  for  output...'  ); 

energies  =  frames . energy; 
intensity  =  f rames .intensity; 


save 

' plott .mat ' 

save 

' plotr .mat ' 

save 

' distances .mat 

save 

' energies .mat ' 

save 

' intensity . mat 

plott; 
plotr; 
distances ; 
energies ; 
intensity; 


display  execution  times 


disp(  'analysis  complete'  ); 
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timeDiff  =  cputime  -  startCPU; 

disp(  sprintf (  '  total  client  cpu  time  =  %9.3f  sec.',  timeDiff  )); 
disp(  sprintf (  '  total  elapsed  time  =  %9.3f  sec.',  toe  )); 


o, 

o 

o, 

o 

o, 

o 


generate  new  plot  frame 


function  frames  =  generateFrame (  figno,  plotting,  field,  frames  ) 


figno  figure  number 

plotting  plotting  control  parameters 

field  field  value  matrix 

frames  updated  stucture 

intensity  =  intensity  matrices  (per  figure) 
energy  =  pulse  energy  (per  figure) 

lastFrame  =  last  frame  number 


frames . lastFrame 


=  figno; 


calculate  pulse  energy 


absField 

sqrField 

powerMatrix 

first Integral 

secondlntegral 

frames .energy (figno) 


abs (  field  ) ; 
absField.* absField; 

sqrField. *plotting .mask; 

trapz (  powerMatrix  ) ; 

trapz (  plotting . radii ,  f irstlntegral '  ); 
secondlntegral *plotting .factor; 


o, _ 

o 

o, 

o 

%  case  A:  wide  view 

o, 

o 


if (  plotting . zoomOption  ==  0  ) 
i4  =  plotting . itplotl ; 

for  iplot  =  1 : plotting . Ntplot 


i4 

=  i4  +  4; 

i3 

=  i4  -  1; 

i2 

=  i3  -  1; 

il 

=  i2  -  1; 

jplot 

=  1; 
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for  j4  =  4 : 4 : plotting . Nradial 
j  3  =  j4  -  1; 

j  2  =  j  3  -  1  ; 

jl  =  j2  -  1; 

magnit  =  sqrField (il,  j 1)  ; 

magnit  =  magnit  +  sqrField ( i2 , j 1 ) ; 
magnit  =  magnit  +  sqrField ( i3 , j 1 ) ; 
magnit  =  magnit  +  sqrField ( i4 , j 1 ) ; 

magnit  =  magnit  +  sqrField (il, j2) ; 
magnit  =  magnit  +  sqrField ( i2 , j 2 ) ; 
magnit  =  magnit  +  sqrField ( i3 , j 2 ) ; 
magnit  =  magnit  +  sqrField ( i4 , j 2 ) ; 

magnit  =  magnit  +  sqrField ( il , j 3 ) ; 
magnit  =  magnit  +  sqrField ( i2 , j 3 ) ; 
magnit  =  magnit  +  sqrField ( i3 , j 3 ) ; 
magnit  =  magnit  +  sqrField ( i4 , j 3 ) ; 

magnit  =  magnit  +  sqrField ( il , j 4 ) ; 
magnit  =  magnit  +  sqrField ( i2 , j 4 ) ; 
magnit  =  magnit  +  sqrField ( i3 , j 4 ) ; 
magnit  =  magnit  +  sqrField(i4,j4); 

frames . intensity (jplot, iplot, figno)  =  0 . 0625*magnit; 
jplot  =  jplot  +  1; 

end 

end 

o, _ 

o 

o, 

o 

%  case  B:  zoomed  view 

o, 

o 

else 

iplot  =  1; 

for  itime  =  plotting . itplotl : plotting . itplot2 
for  jplot  =  1 : plotting . Nrplot 

frames . intensity (jplot, iplot, figno)  =  sqrField (itime, jplot) ; 

end 

iplot  =  iplot  +  1; 

end 

end 


end 


145 


B.2.7  Part  A  of  the  the  “slower  model”  (highResNLSfullPartA) 


For  the  first  part  of  the  “slower  model”,  the  software  uses  highResNLSfullPartA.m: 

function  returned  =  highResNLSfullPartA (  data  ) 

9'*****************'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

part  A  of  propagation  step  --  distributed 
(DHT  processing) 

14  Sep  2006  version  5.5 
for 

MATLAB  7  Distributed  Computing 

input  parameters 

input  =  original  field  values 

returned  data 

o, 

o 

%'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  pre-allocate  arrays 

o, 

o 


o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

%  data 

o, 

o 

o, 

o 

%  returned 


returned . firstDeriv  =  complex ( 1 . 0 , 0 . 0 ).* zeros (  data. size,  data.Nradial  ); 
returned . secondDer iv  =  complex ( 1 . 0 , 0 . 0 ).* zeros (  data. size,  data.Nradial  ); 


o, _ 

o 

o, 

o 

%  radial  derivatives 

o, 

o 


for  item  =  l:data.size 

radialBuf fer2  =  data . input ( item,  ; 

scaled  =  radialBuf fer2 . *data .mlrecip; 

transformed  =  data . XM*scaled; 

o, 

o 

%  first  derivative 

O, 

o 


onceScaled 
f irstTrans 
f irstTerm 

returned . firstDeriv ( item, :) 

O, 

o 


=  data . derivOp . *transf ormed; 

=  data . XM*onceScaled; 

=  f irstTrans . *data .ml ; 

=  (data . divRad . *firstTerm) . ' ; 


146 


o\°  o\°  o\° 


second  derivative 


o, 

o 

twiceScaled  =  data . derivOp . *onceScaled; 

secondTrans  =  data . XM*twiceScaled; 

secondTerm  =  secondTrans . *data .ml ; 

returned . secondDeriv ( item, : )  =  secondTerm. ' ; 

end 

end 
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B.2.8  Part  B  of  the  the  “slower  model”  (highResNLSfullPartB) 


For  the  second  part  of  the  “slower  model”,  the  software  uses  highResNLSfullPartB.m: 

function  returned  =  highResNLSfullPartB (  data  ) 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

part  B  of  propagation  step  --  distributed 
(FFT  processing) 

14  Sep  2006  version  5.5 
for 

MATLAB  7  Distributed  Computing 

input  parameters 

returned  stucture 

adjusted  field  values 

o, 

o 

o, 

o 

<£-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  pre-allocate  arrays 

o, 

o 


o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

%  data 

o, 

o 

%  returned 

o, 

o 


returned  =  complex ( 1 . 0 , 0 . 0 ).* zeros (  data . Ntemporal ,  data. size  ); 


o, 

o 

o, 

o 


o, 

o 

o, 

o 

propagation  step 

for  item  =  l:data. 
timeBuf fer 

.  size 

=  data . input ( : , item) ; 

timeBuf ferSq 

=  timeBuf fer . *timeBuf fer  ; 

con j  Buffer 

=  conj (  timeBuffer  ) ; 

con j  Prod 

=  con j Buff er . *timeBuf fer ; 

realpart 

imagpart 

=  real (  timeBuffer  ) ; 

=  imag (  timeBuffer  ); 

abslntens 

=  (realpart . *realpart)  +  (imagpart 

intensProd 

=  abslntens . *timeBuf fer; 

f irstRadDrv 

secondRadDrv 

=  data . firstRadDeriv (:, item) ; 

=  data . secondRadDeriv (:, item) ; 

*imagpart)  ; 
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forward  FFT 


O, 

o 


transField 

= 

fft  ( 

timeBuffer  )  ; 

transCon j 

= 

fft  ( 

conjBuffer  )  ; 

translntens 

= 

fft  ( 

abslntens  ) ; 

transFirstR 

= 

fft  ( 

firstRadDrv  ) ; 

transSecondR 

= 

fft  ( 

secondRadDrv  ) ; 

O, 

o 

%  apply  derivative  operators 

o, 

o 

f irstDerivFr 

— 

data . 

firstDeriv. *transField; 

secondDerivFr 

= 

data . 

secondDeriv. *transField; 

thirdDerivFr 

— 

data . 

thirdDeriv. *transField; 

con j  DerivFr 

= 

data . 

firstDeriv. *transConj ; 

intensDerivFr 

= 

data . 

firstDeriv. *translntens; 

derivFirstFr 

= 

data . 

firstDeriv. *transFirstR; 

derivSecondFr 

= 

data . 

firstDeriv. *transSecondR 

Q, 

O 

%  reverse  FFT 

O, 

o 

f irstDeriv 

_ 

ifft  ( 

firstDerivFr  ) ; 

secondDeriv 

= 

ifft  ( 

secondDerivFr  )  ; 

thirdDeriv 

= 

ifft  ( 

thirdDerivFr  )  ; 

con j  Deriv 

= 

ifft  ( 

con j DerivFr  ); 

intensDeriv 

= 

ifft  ( 

intensDerivFr  )  ; 

derivFirstRD 

= 

ifft  ( 

derivFirstFr  ) ; 

derivSecondRD 

= 

ifft  ( 

derivSecondFr  ) ; 

o, 

o 

%  updated  field  envelope 

o, 

o 


parti 

part2 

part3 

part4 

part5 

part6 

part7 


=  data . factorl . *secondDeriv; 

=  data . f actor2 . *thirdDeriv; 

=  data . f actor3 . *timeBuf fer ; 

=  data . factor4 . *intensProd; 

=  data .factor5.*conj  Prod . *f irstDeriv; 

=  data . fact or 6 . *timeBuf ferSq . * con j  Deri v 
=  data . fact or 7 . *timeBuf fer . *intensDeriv 


partDl 

partD2 

partD3 

partD4 


=  data . scalel . *f irstRadDrv; 

=  data . scale2 . *secondRadDrv; 

=  data . scale3 . *derivFirstRD; 

=  data . scale4 . *derivSecondRD; 
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deltaQ 


newField 


save  results 


=  parti  -  part2  -  part3 
+  part4  -  part5  -  part6 
+  partDl  +  partD2  -  partD3 

=  timeBuffer  +  deltaQ; 


part7 

partD4 


returned ( : , item) 


newField; 


B.2.9  Pulse  animation  plotting  (highResNLSplotl) 


For  pulse  animation  plotting,  the  software  uses  highResNLSplotl.m: 

function  varargout  =  highResNLSplotl (varargin) 

9'*****************'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  pulse  animation  plotting 

O, 

o 

%  14  Sep  2006  version  5.5 

%  for 

%  MAT LAB  7 


input  files 


plott . mat 
plotr . mat 
distances .mat 
energies .mat 
intensity .mat 


time  axis  grid 
radial  axis  grid 
propagation  distances 
pulse  energies 
intensity  surfaces 


Control  Parameters  (inputs) 


plotting  parameters: 
perspective 

1  for  incoming  perspective 

2  for  outgoing  perspective 

3  for  target  perspective 
animation  option 

0  =  off  (multiple  figures  displayed) 

1  =  on  with  compression  (low-quality  AVI  file  generated) 

(single  figure  updated) 

2  =  on  with  no  compression  (high-quality  AVI  file  generated) 

(single  figure  updated) 


halfview 
0 
1 

export  flag 

1  for  on 
0  for  off 


for  full-surface  view 
for  half-surface  view 


(jpeg  files  generated  per  frame  or  figure) 


main  program 


format  long; 
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%  Begin  initialization  code  -  DO  NOT  EDIT 
gui  Singleton  =  1; 

gui  State  =  struct (  'gui  Name',  mfilename,  ... 

'gui  Singleton',  gui  Singleton,  ... 

'gui  OpeningFcn ' ,  QhighResNLSplotl  OpeningFcn,  ... 

'gui  OutputFcn ' ,  QhighResNLSplotl  OutputFcn,  ... 

' gui_LayoutFcn ' ,  [],  ... 

' gui_Callback '  ,  []  ); 

if (  nargin  &  isstr(  varargin{l}  )  ) 

gui  State. gui  Callback  =  str2func(  varargin{l}  ); 

end 

if (  nargout  ) 

[varargout { 1 : nargout } ]  =  gui  mainfcn (  gui  State,  vararginf : }  ); 

else 

gui  mainfcn (  gui  State,  varargin) : }  ); 

end 

%  End  initialization  code  -  DO  NOT  EDIT 

O _ 

O - 

o, 

o 

%  executes  just  before  highResNISplot  is  made  visible. 

o, 

o 

function  highResNLSplotl  OpeningFcn (  hObject,  eventdata,  handles,  varargin  ) 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  varargin  command  line  arguments  to  highResNISplot  (see  VARARGIN) 

O, 

o 

%  This  function  has  no  output  args,  see  OutputFcn. 

O, 

o 

%  Choose  default  command  line  output  for  highResNISplot 
handles . output  =  hObject; 

%  Update  handles  structure 
guidata (  hObject,  handles  ); 

%  UIWAIT  makes  highResNISplot  wait  for  user  response  (see  UIRESUME) 

%  uiwait (handles . figurel ) ; 

o _ 

O - 

o, 

o 

%  outputs  from  this  function  are  returned  to  the  command  line. 

O, 

o 

function  varargout  =  highResNLSplotl  OutputFcn (  hObject,  eventdata,  handles  ) 
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%  hObject 
%  eventdata 
%  handles 

o, 

o 

%  varargout 


handle  to  figure 

reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
structure  with  handles  and  user  data  (see  GUIDATA) 

cell  array  for  returning  output  args  (see  VARARGOUT); 


%  Get  default  command  line  output  from  handles  structure 


varargout! 1}  =  handles . output; 


incoming  perspective 


function  incoming  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  incoming  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value')  returns  toggle  state  of  incoming 

disp(  'incoming  perspective  selected'  ); 

set (  handles . target,  'Value',  0  ); 
set (  handles . outgoing,  'Value',  0  ); 


o, 

o 

o, 

o 

o, 

o 


outgoing  perspective 


function  outgoing_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  outgoing  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value ' )  returns  toggle  state  of  outgoing 

disp(  'outgoing  perspective  selected'  ); 

set (  handles . target,  'Value',  0  ); 
set (  handles . incoming,  'Value',  0  ); 


o, 

o 

o, 

o 

o, 

o 


target  perspective 


function  target_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  target  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value ' )  returns  toggle  state  of  target 
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disp(  'target  perspective  selected'  ); 

set (  handles . incoming,  'Value',  0  ); 
set (  handles . outgoing,  'Value',  0  ); 


o, 

o 

o, 

o 

o, 

o 


uncompressed  animation 


function  uncompressed_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  uncompressed  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value')  returns  toggle  state  of  uncompressed 

disp(  'uncompressed  selected'  ); 

set (  handles . compressed,  'Value',  0  ); 
set (  handles . noanimation,  'Value',  0  ); 


o, 

o 

o, 

o 


o, 

o 


compressed  animation 


function  compressed  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  compressed  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObj ect, ' Value ' )  returns  toggle  state  of  compressed 

disp(  'compressed  selected'  ); 

set (  handles . uncompressed,  'Value',  0  ); 
set (  handles . noanimation,  'Value',  0  ); 


o, 

o 

o, 

o 

o, 

o 


no  animation  (multiple  figures) 


function  noanimation  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  noanimation  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value')  returns  toggle  state  of  noanimation 

disp(  'noanimation  selected'  ); 

set (  handles . uncompressed,  'Value',  0  ); 
set (  handles . compressed,  'Value',  0  ); 
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export  toggle 


function  export  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  export  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value ' )  returns  toggle  state  of  export 


disp(  'export  selected'  ); 


o _ 

O - 

%  half-view  toggle 

o _ 

O - 


function  halfview_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  halfview  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value')  returns  toggle  state  of  halfview 


disp(  'halfview  selected'  ); 


o _ 

O - 

%  respond  to  button  press  (activate  plotting) 

O _ 

O - 


function  plot  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  plot  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 


set  ( 

handles .closeall. 

' BackgroundColor '  , 

'  red ' 

)  ; 

set  ( 

handles .closeall. 

' Enable '  , 

'off' 

)  ; 

set  ( 

handles .plot. 

' BackgroundColor ' , 

'  red ' 

)  ; 

set  ( 

handles .plot. 

' Enable ' , 

'off' 

) ; 

disp(  'plot  activated'  ); 
format  long; 

O, _ 

o 

o, 

o 

%  parameters  for  plotting 

O, 

o 

if (  get (  handles . target,  'Value'  )  >  0  ) 
perspective  =  3; 

elseif  (  get  (  handles . outgoing,  'Value'  )  >  0  ) 
perspective  =  2; 
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else 


perspective  =  1; 

end 

if (  get (  handles . uncompressed,  'Value'  )  >  0  ) 
animation  =  2; 

elseif (  get (  handles . compressed,  'Value'  )  >  0  ) 


animation 

=  1; 

animation 

=  0; 

end 

halfview  =  get (  handles . halfview,  'Value'  ); 
exportFlag  =  get (  handles . export,  'Value'  ); 


%d',  perspective 
%d',  animation 
%d ' ,  halfview 
%d',  exportFlag 


disp(  'plotting  parameters:'  ); 

disp(  sprintf (  '  perspective  ( l=incoming, 2=outgoing, 3=target)  = 

) )  ; 

disp(  sprintf (  '  animation  (2=uncompr, l=compr, 0=of f )  = 

) )  ; 

disp(  sprintf (  '  halfview  flag  ( l=half , 0=full )  = 

) )  ; 

disp(  sprintf (  '  export  flag  (l=on,0=off)  = 


o, 

o 

o, 

o 


o, 

o 


load  data 


disp  ( 

' loading  output 

.  data . .  .  ' 

load 

' plott .mat ' 

plott; 

load 

' plotr . mat ' 

plotr; 

load 

' distances .mat ' 

distances ; 

load 

' energies . mat ' 

energies ; 

load 

' intensity . mat ' 

intensity; 

o, 

o 

%  radial  axis 

o, 

o 


lengths  =  size(  plotr  ); 

Nrdata  =  lengths (1); 


radialMax  =  plotr (Nrdata) ; 

if (  plotr(l)  ==  0.0  ) 

Nrplot  =  2*Nrdata  -  1; 
else 

Nrplot  =  2*Nrdata; 

end 


O, 

o 

%  temporal  axis 

O, 

o 
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lengths  =  size(  plott  ); 

Ntplot  =  lengths (1); 

timeMax  =  plott (Ntplot) ; 


o, 

o 

%  adjust  plotting  units 

O, 

o 


if (  radialMax  <  (10A(-3))  ) 

plotr  =  plotr* (10A6) ; 

radialMax  =  radialMax* (10A6) ; 

radialLabel  =  'radius  (micron) '; 

elseif  (  radialMax  <1.0  ) 

plotr  =  plotr* (10A3) ; 

radialMax  =  radialMax* ( 1 0  A  3 ) ; 

radialLabel  =  'radius  (mm) '; 

else 

radialLabel  =  'radius  (m) '; 

end 

if (  timeMax  <  (10A (—12) )  ) 

plott  =  plott* (10A15) ; 

timeMax  =  timeMax* ( 1 0  A 1 5 ) ; 
timeLabel  =  'time  (fs) '; 

elseif (  timeMax  <  (10A(-9))  ) 

plott  =  plott* (10A12) ; 

timeMax  =  timeMax* (10 A 12) ; 
timeLabel  =  'time  (ps) '; 

elseif (  timeMax  <  (10A(-6))  ) 

plott  =  plott* (10A9) ; 

timeMax  =  timeMax* (10 A 9)  ; 
timeLabel  =  'time  (ns) '; 

else 

timeLabel  =  'time  (sec) '; 

end 


if ( (  Ntplot  ==  128  ) 

I  1 

(  Ntplot  == 

256  )  ) 

verticalLabel  = 

'  rel . 

intensity 

(4x4  Avg 

else 

verticalLabel  = 

'  rel . 

intensity ' 

r 

end 

o, 

o 

%  mirror  r  about  z 

axis 

o, 

o 


lengths  =  size(  intensity  ); 
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Numplots  =  lengths  (3); 

lastFigure  =  0; 

for  figno  =  1: Numplots 

if (  energies (figno)  ==  0.0  ) 
break; 

end 

lastFigure  =  figno; 

end 

Numplots  =  lastFigure; 

fullr  =  zeros (  Nrplot,  1  ) ; 

fullSurf  =  zeros (  Nrplot,  Ntplot,  Numplots  ) ; 

for  figno  =  1: Numplots 
for  iplot  =  l:Ntplot 
jplot  =  1; 

for  jradial  =  Nrdata:-l:l 

fullr(jplot)  =  -plotr ( j radial ) ; 

fullSurf (jplot, iplot, figno)  =  intensity (jradial, iplot, figno) ; 

jplot  =  jplot  +  1; 

end 

if (  plotr(l)  ==  0.0  ) 

for  jradial  =  2:Nrdata 
fullr (jplot) 

fullSurf (jplot, iplot, figno) 
jplot 

end 
else 

for  jradial  =  l:Nrdata 
fullr (jplot) 

fullSurf (jplot, iplot, figno) 
jplot 

end 

end 

end 

end 

set (  handles . closeall ,  'UserData',  Numplots  ); 

o, _ 

o 

%  output  initial  pulse 

O, 


=  plotr ( j radial )  ; 

=  intensity (j  radial, iplot,  figno)  ; 
=  jplot  +  1; 

=  plotr ( j radial ) ; 

=  intensity (j  radial, iplot, figno) ; 
=  jplot  +  1; 
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disp(  'plotting...'  ) 


propagDist  =  0.0; 
figno  =  1; 

handle  =  figure (  1  ) ; 

pulseEnergy  =  energies (figno) ; 

etext  =  sprintf  (  'pulse  energy  (FWHM)  = 

text  =  sprintf (  'prop,  distance  =  0.0, 

axis (  [-timeMax  timeMax  -radialMax  radialMax] 


%8.6f  J',  pulseEnergy  ) 
%s ' ,  etext  ) ; 

)  ; 


title (  text,  ' FontWeight '  ,  'bold'  ); 

xlabel (  timeLabel,  'FontWeight',  'bold'  ); 

ylabel (  radialLabel,  'FontWeight',  'bold'  ); 

zlabel (  verticalLabel ,  'FontWeight',  'bold'  ); 

if (  perspective  ==  3  ) 
view (  [90  0]  )  ; 
elseif (  perspective  ==  2  ) 
view (  [ -50  30 ]  )  ; 
else 

view (  [50  30 ]  )  ; 

end 


pause  (  0.1  )  ; 

hold  on; 
grid  on; 

if (  halfview  ==  1  ) 

surfl (  plott,  plotr,  intensity figno)  ); 
else 

surfl (  plott,  fullr,  fullSurf figno)  ); 

end 

shading  interp; 
colormap (  gray  ) ; 
hold  of f ; 

refresh (  1  ) ; 

if (  animation  ~=  0  ) 

frames (figno)  =  getframe (  handle  ); 

end 

if (  exportFlag  ~=  0  ) 

filename  =  sprintf (  'figure%d.jpg',  figno  ); 
print (  handle,  '-djpeg',  filename  ); 

end 
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o\°  o\°  o\° 


propagation  loop 


for  figno  =  2:Numplots 

propagDist  =  distances (figno)  ; 

text  =  sprintf (  '  frame  %d  of  %d',  figno,  Numplots  ); 

set  (  handles . status ,  'String',  text  ); 

if (  propagDist  >  0  ) 

pulseEnergy  =  energies (figno) ; 

if (  animation  ==  0  ) 

handle  =  figure (  figno  ) ; 

else 

elf  ; 

end 

if (  propagDist  <  (10A-6)  ) 

dtext  =  sprintf (  'prop,  distance  =  %5.3f  nm',  100000000 . 0*propagDist  ) 
elseif (  propagDist  <  (10A-3)  ) 

dtext  =  sprintf  (  'prop,  distance  =  %5.3f  micron',  1000000 . 0*propagDist 

)  ; 

elseif (  propagDist  <1.0  ) 

dtext  =  sprintf (  'prop,  distance  =  %5.3f  mm',  1000 . 0*propagDist  ); 
else 

dtext  =  sprintf (  'prop,  distance  =  %5.3f  m',  propagDist  ); 

end 

etext  =  sprintf (  'pulse  energy  (FWHM)  =  %8.6f  J',  pulseEnergy  ); 

text  =  sprintf (  '%s,  %s',  dtext,  etext  ); 

axis (  [-timeMax  timeMax  -radialMax  radialMax]  ); 

title (  text,  ' FontWeight ' ,  'bold'  ); 

xlabel (  timeLabel,  'FontWeight',  'bold'  ); 

ylabel (  radialLabel,  'FontWeight',  'bold'  ); 

zlabel (  verticalLabel ,  'FontWeight',  'bold'  ); 

if (  perspective  ==  3  ) 


view  ( 

[90  0]  )  ; 

elseif ( 

perspective 

view  ( 

[-50  30]  ) 

else 

view  ( 

[50  30]  ) ; 

end 

hold  on; 

grid  on; 
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o\°  o\°  o\° 


if (  halfview  ==  1  ) 

surfl (  plott,  plotr,  intensity figno)  ); 
else 

surfl (  plott,  fullr,  fullSurf figno)  ); 

end 

shading  interp; 
colormap (  gray  ) ; 
hold  off; 

if (  animation  ==  0  ) 
refresh (  figno  ); 

else 

drawnow; 
refresh (  1  ) ; 

frames (figno)  =  getframe (  handle  ); 

end 

if (  exportFlag  ~=  0  ) 

filename  =  sprintf (  'figure%d.jpg',  figno  ); 
print (  handle,  ' -djpeg',  filename  ); 

end 

end 

end 


save  movie 


if (  animation  ~=  0  ) 

disp(  'saving  to  avi . . . '  ); 

text  =  sprintf (  'saving  to  AVI  (%d  frames) ',  Numplots  ); 

set (  handles . status ,  'String',  text  ); 

if (  perspective  ==  3  ) 
if (  animation  ==  2  ) 


movie2avi ( 
else 

frames , 

' target . avi ' , 

' FPS  '  , 

3, 

'COMPRESSION' , 

' None ' 

movie2avi ( 
QUALITY',  100  ); 

frames , 

' target . avi ' , 

' FPS  '  , 

3, 

'COMPRESSION' , 

' Indeo5 

end 

disp(  '  target.avi  saved'  ); 

set (  handles . status ,  'String',  'target.avi  saved'  ); 

elseif (  perspective  ==  2  ) 
if (  animation  ==  2  ) 

movie2avi (  frames,  'outgoing.avi',  ' FPS ' ,  3,  'COMPRESSION',  'None'  ) 
else 


161 


movie2avi (  frames,  'outgoing.avi',  ' FPS ' ,  3,  'COMPRESSION',  ' Indeo5 ' 
'QUALITY',  100  ) ; 
end 

disp(  '  outgoing.avi  saved'  ); 

set (  handles . status ,  'String',  'outgoing.avi  saved'  ); 
else 

if (  animation  ==  2  ) 


movie2avi ( 
else 

frames , 

' incoming . avi ' , 

' FPS  '  , 

3, 

'COMPRESSION' , 

' None ' 

movie2avi ( 
QUALITY',  100  ); 

frames , 

'  incoming . avi ' , 

' FPS  '  , 

3, 

'COMPRESSION' , 

' Indeo5 

end 

disp(  '  incoming.avi  saved'  ); 

set (  handles . status ,  'String',  'incoming.avi  saved'  ); 

end 

end 

set (  handles .plot, 
set (  handles .plot, 
set (  handles . closeall , 
set (  handles . closeall , 

o _ 

o - 

%  respond  to  button  press  (close  all) 

o _ 

O - 


'  BackgroundColor ' ,  'green'  ); 
'  Enable ' ,  ' on '  )  ; 
'BackgroundColor',  'green'  ); 
'  Enable ' ,  ' on '  )  ; 


function  closeall_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  closeall  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

if (  get (  handles . noanimation,  'Value'  )  >  0  ) 

%  close  all; 

Numplots  =  get (  handles . closeall ,  'UserData'  ); 

for  figno  =  lrNumplots 

handle  =  figure  (  figno  ) ; 

close  (  handle  ) ; 

end 

else 

handle  =  figure  (  1  ) ; 
close  (  handle  ) ; 

end 

set (  handles . closeall ,  'Enable',  'off'  ); 
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B.2.10  Composite  profile  plotting  (highResNLSplot2) 


For  the  composite  profile  displays,  the  software  uses  highResNLSplot2.m: 


function  varargout  =  highResNLSplot2 (varargin) 

9'****************'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  cummulative  profile  plotting 

O, 

o 

%  14  Sep  2006  version  5.5 

%  for 

%  MAT LAB  7 


plott . mat 
plotr .mat 
distances .mat 
energies .mat 
intensity .mat 


input  files 


time  axis  grid 
radial  axis  grid 
propagation  distances 
pulse  energies 
intensity  surfaces 


o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 

o, 

o 


Control  Parameters  (inputs) 


plotting  parameters: 
perspective 

1  for  incoming  perspective 

2  for  outgoing  perspective 

3  for  target  perspective 


O, 

o 

o, 

o 

o, 

o 


main  program 


format  long; 


%  Begin  initialization  code  -  DO  NOT  EDIT 


gui  Singleton  =  1; 

mf ilename, 
gui  Singleton, 

@highResNLSplot2  OpeningFcn, 
@highResNLSplot2  OutputFcn, 
[], 

[] 

if (  nargin  &  isstr (varargin { 1 } )  ) 


gui_State  = 


struct (  ' gui 
'  gui 
'  gui 
'  gui 
'  gui 
'  gui 


Name ' , 
Singleton ' , 
OpeningFcn ' , 
OutputFcn ' , 
LayoutFcn ' , 
Callback ' , 
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gui  State . gui^Callback  =  str2func(  varargin{l}  ); 

end 

if (  nargout  ) 

[varargout { 1 : nargout } ]  =  gui  mainfcn (  gui  State,  varargin{ : }  ); 

else 

gui  mainfcn (  gui  State,  varargin) : }  ); 

end 

%  End  initialization  code  -  DO  NOT  EDIT 

O _ 

O - 

o, 

o 

%  executes  just  before  highResNLSplot2  is  made  visible. 

o, 

o 

function  highResNLSplot2  OpeningFcn (  hObject,  eventdata,  handles,  varargin  ) 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  varargin  command  line  arguments  to  highResNISplot  (see  VARARGIN) 

O, 

o 

%  This  function  has  no  output  args,  see  OutputFcn . 

O, 

o 

%  Choose  default  command  line  output  for  highResNISplot 
handles . output  =  hObject; 

%  Update  handles  structure 
guidata (  hObject,  handles  ); 

%  UIWAIT  makes  highResNISplot  wait  for  user  response  (see  UIRESUME) 

%  uiwait (handles . figurel )  ; 

plot_Callback (  handles  ); 

o _ 

O - 

o, 

o 

%  outputs  from  this  function  are  returned  to  the  command  line. 

O, 

o 

function  varargout  =  highResNLSplot2  OutputFcn  (  hObject,  eventdata,  handles  ) 
%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

O, 

o 

%  varargout  cell  array  for  returning  output  args  (see  VARARGOUT); 

%  Get  default  command  line  output  from  handles  structure 
varargout{l}  =  handles .output; 
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incoming  perspective 


function  incoming_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  incoming  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value')  returns  toggle  state  of  incoming 

disp(  'incoming  perspective  selected'  ); 

set (  handles . target,  'Value',  0  ); 
set (  handles . outgoing,  'Value',  0  ); 

plot_Callback (  handles  ); 


outgoing  perspective 


function  outgoing_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  outgoing  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value ' )  returns  toggle  state  of  outgoing 

disp(  'outgoing  perspective  selected'  ); 

set (  handles . target,  'Value',  0  ); 
set (  handles . incoming,  'Value',  0  ); 

plot_Callback (  handles  ); 


o, 

o 


target  perspective 


function  target  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  target  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value ' )  returns  toggle  state  of  target 

disp(  'target  perspective  selected'  ); 

set (  handles . incoming,  'Value',  0  ); 
set (  handles . outgoing,  'Value',  0  ); 

plot_Callback (  handles  ); 
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maximum  (time)  value  mode 


function  maxVal  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  maxVal  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'maximum  value  mode  selected'  ); 


set (  handles . integrated,  'Value',  0  ); 


plot_Callback (  handles  ); 


o _ 

O - 

%  integrated  mode 

o _ 

O - 


function  integrated  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  integrated  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'integrated  mode  selected'  ); 


set (  handles . maxVal ,  'Value',  0  ); 


plot_Callback (  handles  ); 


o _ 

O - 

%  respond  to  button  press  (activate  plotting) 

O _ 

O - 


function  plot  Callback (  handles  ) 

%  hObject  handle  to  plot  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'plot  activated'  ); 

format  long; 


O, _ 

o 

o, 

o 

%  parameters  for  plotting 

O, 

o 


if (  get (  handles . target,  'Value'  )  >  0  ) 
perspective  =  3; 

elseif  (  get  (  handles . outgoing,  'Value'  )  >  0  ) 
perspective  =  2; 
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else 


perspective  =  1; 

end 

if (  get (  handles . integrated,  'Value'  )  >  0  ) 
integrated  =  1; 

else 

integrated  =  0; 

end 

disp(  'plotting  parameters:'  ); 

disp(  sprintf (  '  perspective  =  %d',  perspective  )  ); 
disp(  sprintf (  '  integration  =  %d',  integrated  )  ); 


O, 

o 

o, 

o 


o, 

o 


load  data 


disp  ( 

' loading  output 

.  data . .  .  ' 

load 

' plott . mat ' 

plott; 

load 

' plotr .mat ' 

plotr; 

load 

' distances . mat ' 

distances ; 

load 

' energies . mat ' 

energies ; 

load 

' intensity .mat ' 

intensity; 

o, 

o 

%  radial  axis 

o, 

o 


lengths  =  size(  plotr  ); 

Nrdata  =  lengths (1); 

radialMax  =  plotr (Nrdata) ; 

if (  plotr(l)  ==  0.0  ) 

Nrplot  =  2*Nrdata  -  1; 
else 

Nrplot  =  2*Nrdata; 

end 

O, 

o 

%  distance  axis 

o, 

o 


lengths  =  size(  distances  ); 

Nzplot  =  lengths (1); 

lastFigure  =  0; 

for  figno  =  l:Nzplot 

if (  energies (figno)  ==  0.0  ) 
break; 

end 
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lastFigure  =  figno; 

end 

distMax  =  distances (lastFigure) ; 

o, 

o 

%  temporal  axis 

O, 

o 

lengths  =  size(  plott  ); 

Ntplot  =  lengths (1); 

O, 

o 

%  adjust  plotting  units 

O, 

o 

if (  radialMax  <  (10A(-3))  ) 

plotr  =  plotr* (10A6) ; 

radialMax  =  radialMax* (10A6) ; 
radialLabel  =  'radius  (micron) '; 

elseif (  radialMax  <1.0  ) 

plotr  =  plotr*  (10A3) ; 

radialMax  =  radialMax* ( 1 0  A  3 ) ; 
radialLabel  =  'radius  (mm) '; 

else 

radialLabel  =  'radius  (m) '; 

end 

distLabel  =  'distance  (m) '; 

O, 

o 

%  mirror  r  about  z  axis 

O, 

o 

lengths  =  size(  intensity  ); 

Numplots  =  lengths (3); 

lastFigure  =  0; 

for  figno  =  l:Numplots 

if (  energies (figno)  ==  0.0  ) 
break; 

end 

lastFigure  =  figno; 

end 

fullr  =  zeros  (  Nrplot,  1  )  ; 
jplot  =  1; 

for  jradial  =  Nrdata:-l:l 
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fullr(jplot)  =  -plotr ( j radial )  ; 
jplot  =  jplot  +  1; 

end 

if (  plotr(l)  ==  0.0  ) 

for  j radial  =  2:Nrdata 

fullr(jplot)  =  plotr ( j radial ) ; 

jplot  =  jplot  +  1; 

end 

else 

for  j radial  =  l:Nrdata 

fullr(jplot)  =  plotr ( j radial )  ; 

jplot  =  jplot  +  1; 

end 

end 

fullSurf  =  NaN (  Nrplot,  Nzplot  ); 

if (  integrated  >  0  ) 
maxval  =  0.0; 

values  =  zeros  (  Ntplot,  1  )  ; 

deltat  =  plott(2)  -  plott(l); 

for  figno  =  l:lastFigure 
jplot  =  1; 

for  jradial  =  Nrdata:-l:l 

values  =  intensity (j radial, :, figno) ; 

fullSurf (jplot, figno)  =  deltat*trapz (  values  ); 

if (  fullSurf (jplot, figno)  >  maxval  ) 
maxval  =  fullSurf (jplot,  figno)  ; 

end 

jplot  =  jplot  +  1; 

end 

if (  plotr(l)  ==  0.0  ) 

for  jradial  =  2:Nrdata 

values  =  intensity (j radial, :, figno) 

fullSurf (jplot, figno)  =  deltat*trapz (  values  ); 

jplot  =  jplot  +  1; 

end 

else 

for  jradial  =  l:Nrdata 

values  =  intensity (j radial, :, figno) 
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fullSurf ( jplot, figno)  =  deltat*trapz (  values  ); 
jplot  =  jplot  +  1; 

end 

end 

end 

normalizer  =  1.0/maxval; 

for  figno  =  l:lastFigure 
for  jradial  =  l:Nrplot 

fullSurf (jradial, figno)  =  fullSurf (j radial, figno) *normalizer 

end 

end 

else 

for  figno  =  l:lastFigure 
jplot  =  1; 

for  jradial  =  Nrdata:-l:l 
plotted  =  0.0; 

for  iplot  =  l:Ntplot 

value  =  intensity (jradial, iplot, figno) ; 

if (  value  >  plotted  ) 
plotted  =  value; 

end 

end 

fullSurf (jplot, figno)  =  plotted; 
jplot  =  jplot  +  1; 

end 

if (  plotr(l)  ==  0.0  ) 

for  jradial  =  2:Nrdata 
plotted  =  0.0; 

for  iplot  =  l:Ntplot 

value  =  intensity (jradial, iplot,  figno) ; 

if (  value  >  plotted  ) 
plotted  =  value; 

end 

end 

fullSurf (jplot, figno)  =  plotted; 
jplot  =  jplot  +  1; 

end 

else 

for  jradial  =  l:Nrdata 
plotted  =  0.0; 
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for  iplot  =  l:Ntplot 

value  =  intensity (j radial, iplot, figno) ; 

if (  value  >  plotted  ) 
plotted  =  value; 

end 

end 

fullSurf (jplot, figno)  =  plotted; 
jplot  =  jplot  +  1; 

end 

end 

end 

end 

o, _ 

o 

%  output  propagation  silhouette 

O, _ 

o 

disp(  'plotting...'  ); 
figure ( 1 )  ; 
elf  ; 

axis (  [0.0  distMax  -radialMax  radialMax]  ); 

if (  integrated  >  0  ) 

title  (  'normalized  (time)  integrated  propagated  pulse',  ' FontWeight ' ,  'bold' 

)  ; 

if ( (  Ntplot  ==  128  )  I  I  (  Ntplot  ==  256  )) 

verticalLabel  =  'integrated  rel.  intensity  (4x4  Avg.) '; 
else 

verticalLabel  =  'integrated  rel.  intensity'; 

end 

else 

title  (  'maxima  of  propagated  pulse',  'FontWeight',  'bold'  ); 

if ( (  Ntplot  ==  128  )  I  I  (  Ntplot  ==  256  )) 

verticalLabel  =  'maximum  rel.  intensity  (4x4  Avg.) '; 
else 

verticalLabel  =  'maximum  rel.  intensity'; 

end 

end 

xlabel (  distLabel,  'FontWeight',  'bold'  ); 

ylabel (  radialLabel,  'FontWeight',  'bold'  ); 

zlabel (  verticalLabel,  'FontWeight',  'bold'  ); 

if (  perspective  ==  3  ) 
view (  [90  0]  )  ; 
elseif (  perspective  ==  2  ) 
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view (  [ -50  30 ]  )  ; 
else 

view (  [50  30 ]  )  ; 

end 

pause  (  0.1  )  ; 

hold  on; 
grid  on; 

surfl (  distances,  fullr,  fullSurf  ) 

shading  interp; 
colormap (  gray  ) ; 
hold  off; 


refresh (  1  )  ; 


B.2.11  Target-plane  animation  plotting  (highResNLSplot3) 


For  animations  involving  target  plane  energy  patterns,  the  software  uses  highResNLSplot3.m: 

function  varargout  =  highResNLSplot3 (varargin) 

9'******************'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

o, 

o 

%  Laser  Beam  Propagation  through  Air 

O, 

o 

%  target  plane  animation 

O, 

o 

%  14  Sep  2006  version  5.5 

%  for 

%  MAT LAB  7 


input  files 


plott . mat 
plotr . mat 
distances .mat 
energies .mat 
intensity .mat 


time  axis  grid 
radial  axis  grid 
propagation  distances 
pulse  energies 
intensity  surfaces 


Control  Parameters  (inputs) 

plotting  parameters: 
mode 

1  for  maximum  values 

2  for  integrated  values 
animation  option 

0  =  off  (multiple  figures  displayed) 

1  =  on  with  compression  (low-quality  AVI  file  generated) 

(single  figure  updated) 

2  =  on  with  no  compression  (high-quality  AVI  file  generated) 

(single  figure  updated) 


main  program 


format  long; 

%  Begin  initialization  code  -  DO  NOT  EDIT 
gui  Singleton  =  1; 

gui  State  =  struct (  'gui  Name',  mfilename,  ... 

'gui  Singleton',  gui  Singleton,  ... 

'gui  OpeningFcn ' ,  @highResNLSplot3  OpeningFcn,  ... 
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' gui  OutputFcn ' ,  @highResNLSplot3  OutputFcn,  ... 

' gui_LayoutFcn ' ,  [],  ... 

' gui_Callback ' ,  []  ); 

if (  nargin  &  isstr (varargin { 1 } )  ) 

gui  State . gui_Callback  =  str2func(  varargin{l}  ); 

end 

if (  nargout  ) 

[varargout { 1 : nargout } ]  =  gui  mainfcn (  gui  State,  vararginf : }  ); 
else 

gui  mainfcn (  gui  State,  varargin) : }  ); 

end 

%  End  initialization  code  -  DO  NOT  EDIT 

O _ 

O - 

o, 

o 

%  executes  just  before  highResNISplot  is  made  visible. 

o, 

o 

function  highResNLSplot3  OpeningFcn (  hObject,  eventdata,  handles,  varargin  ) 

%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  varargin  command  line  arguments  to  highResNISplot  (see  VARARGIN) 

O, 

o 

%  This  function  has  no  output  args,  see  OutputFcn. 

O, 

o 

%  Choose  default  command  line  output  for  highResNISplot 
handles . output  =  hObject; 

%  Update  handles  structure 
guidata (  hObject,  handles  ); 

%  UIWAIT  makes  highResNISplot  wait  for  user  response  (see  UIRESUME) 

%  uiwait (handles . figurel ) ; 

o _ 

O - 

o, 

o 

%  outputs  from  this  function  are  returned  to  the  command  line. 

O, 

o 

function  varargout  =  highResNLSplot3  OutputFcn (  hObject,  eventdata,  handles  ) 
%  hObject  handle  to  figure 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

O, 

o 

%  varargout  cell  array  for  returning  output  args  (see  VARARGOUT); 

%  Get  default  command  line  output  from  handles  structure 
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varargout{l}  =  handles . output; 


O, 

o 


maximum  (time)  value  mode 


function  maxVal  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  maxVal  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'maximum  value  mode  selected'  ); 


set (  handles . integrated,  'Value',  0  ); 


o _ 

O - 

%  integrated  mode 

o _ 

O - 


function  integrated  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  integrated  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

disp(  'integrated  mode  selected'  ); 


set (  handles . maxVal ,  'Value',  0  ); 


o _ 

O - 

%  uncompressed  animation 

O _ 

O - 


function  uncompressed_Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  uncompressed  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint;  get (hObject, 'Value')  returns  toggle  state  of  uncompressed 

disp(  'uncompressed  selected'  ); 

set (  handles . compressed,  'Value',  0  ); 
set (  handles . noanimation,  'Value',  0  ); 


o, 

o 


compressed  animation 


function  compressed  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  compressed  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
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handles 


structure  with  handles  and  user  data  (see  GUIDATA) 


%  Hint:  get (hObject, 'Value ' )  returns  toggle  state  of  compressed 

disp(  'compressed  selected'  ); 

set (  handles . uncompressed,  'Value',  0  ); 
set (  handles . noanimation,  'Value',  0  ); 

o _ 

O - 

%  no  animation  (multiple  figures) 

O _ 

O - 


function  noanimation  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  noanimation  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

%  Hint:  get (hObject, 'Value')  returns  toggle  state  of  noanimation 

disp(  'noanimation  selected'  ); 

set (  handles . uncompressed,  'Value',  0  ); 
set (  handles . compressed,  'Value',  0  ); 

o _ 

O - 

%  draw  in  2-D  view 

O _ 

O - 


function  draw2D  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  draw2D  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

set (  handles . draw3D,  'Value',  0  ); 

o _ 

O - 

%  draw  in  3-D  view 

O _ 

O - 


function  draw3D  Callback (  hObject,  eventdata,  handles  ) 

%  hObject  handle  to  draw3D  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 

set (  handles . draw2D,  'Value',  0  ); 

o _ 

O - 

%  respond  to  button  press  (activate  plotting) 

O _ 

O - 


function  plot  Callback (  hObject,  eventdata,  handles  ) 
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%  hObject  handle  to  plot  (see  GCBO) 

%  eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
%  handles  structure  with  handles  and  user  data  (see  GUIDATA) 


set  ( 

handles .closeall. 

' BackgroundColor '  , 

'  red ' 

)  ; 

set  ( 

handles .closeall. 

' Enable ' , 

'off' 

)  ; 

set  ( 

handles .plot. 

' BackgroundColor '  , 

'  red ' 

)  ; 

set  ( 

handles .plot. 

' Enable ' , 

'off' 

) ; 

disp(  'plot  activated'  ); 
format  long; 

O, _ 

o 

o, 

o 

%  parameters  for  plotting 

O, 

o 

if (  get (  handles . integrated,  'Value'  )  >  0  ) 
integrated  =  1; 

else 

integrated  =  0; 

end 

if (  get (  handles . uncompressed,  'Value'  )  >  0  ) 
animation  =  2; 

elseif (  get (  handles . compressed,  'Value'  )  >  0  ) 
animation  =  1; 

else 

animation  =  0; 

end 

disp(  'plotting  parameters:'  ); 

disp(  sprintf (  '  integration  =  %d',  integrated  )  )  ; 
disp(  sprintf (  '  animation  =  %d',  animation  )  ); 


O, 

o 

o, 

o 


o, 

o 


load  data 


disp  ( 

' loading  output 

.  data . . .  ' 

load 

'  plott . mat ' 

plott; 

load 

'  plotr .mat ' 

plotr; 

load 

'  distances .mat ' 

distances ; 

load 

'  energies .mat ' 

energies ; 

load 

'  intensity .mat ' 

intensity; 

o, 

o 

%  theta  axis 

o, 

o 


Ntheta  =  361; 


o, 

o 

%  radial  axis 
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o, 

o 


lengths 

=  size (  plotr  ) ; 

Nrdata 

=  lengths ( 1 ) ; 

Nrhalf 

=  f loor (  Nrdata/3  ); 

radialMax 

=  plotr (Nrhalf ) ; 

Nrplot 

=  2*Nrhalf; 

Nraxial 

=  Nrhalf  +  1; 

o, 

o 

%  distance  axis 

o, 

o 

lengths  =  size(  distances  ); 

Nzplot  =  lengths (1); 

O, 

o 

%  temporal  axis 

O, 

o 

lengths  =  size(  plott  ); 

Ntplot  =  lengths (1); 

O, 

o 

%  adjust  plotting  units 

O, 

o 

if (  radialMax  <  (10A(-3))  ) 

plotr  =  plotr* (10A6) ; 

radialMax  =  radialMax* (10A6) ; 
radialLabel  =  'radius  (micron) 

elseif (  radialMax  <1.0  ) 

plotr  =  plotr*  (10A3) ; 

radialMax  =  radialMax* ( 1 0  A  3 ) ; 
radialLabel  =  'radius  (mm) '; 

else 

radialLabel  =  'radius  (m) '; 

end 

distLabel  =  'distance  (m) '; 

O, _ 

o 

o, 

o 

%  find  last  figure 

O, 

o 


lengths  =  size(  intensity  ); 
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Numplots  =  lengths  (3); 

lastFigure  =  0; 

for  figno  =  1: Numplots 

if (  energies (figno)  ==  0.0  ) 

end 

lastFigure  =  figno; 

end 

set (  handles . closeall ,  'UserData',  Numplots  ) 

o, _ 

o 

o, 

o 

%  radial  axis  (mirrored  about  z  axis) 

O, 

o 


fullr  =  zeros  (  Nrplot,  1  )  ; 
jplot  =  1; 

for  jradial  =  Nrhalf:-l:l 

fullr(jplot)  =  -plotr (j radial) ; 

jplot  =  jplot  +  1; 

end 

for  jradial  =  lrNrhalf 

fullr (jplot)  =  plotr ( j radial ) ; 

jplot  =  jplot  +  1; 

end 

o, _ 

o 

o, 

o 

%  surface  of  revolution 

O, 

o 


spun  =  NaN (  Ntheta,  Nraxial,  Nzplot  ); 

for  itheta  =  l:Ntheta 
xx(itheta,l)  =  0.0; 
yy (itheta, 1)  =  0.0; 

end 

jl  ml; 

for  j2  =  2:Nraxial 

rvalue  =  plotr (jl); 

for  itheta  =  lrNtheta 

theta  =  itheta*pi/180 . 0; 
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o\°  o\°  o\° 


xx ( itheta, j 2 )  =  rvalue*cos (  theta  ); 
yy ( itheta, j 2 )  =  rvalue*sin(  theta  ); 

end 

jl  =  j2; 

end 


case  A:  integrated 


if (  integrated  >  0  ) 

zMax  =  1.0; 

verticalLabel  =  'norm,  integrated  intensity'; 

maxval  =  0.0; 

values  =  zeros (  Ntplot,  1  ) ; 

deltat  =  plott(2)  -  plott(l); 

for  figno  =  l:lastFigure 

ji  =  i; 

for  j2  =  2:Nraxial 

values  =  intensity (jl,  figno)  ; 

axval  =  deltat*trapz (  values  ) ; 

if (  axval  >  maxval  ) 
maxval  =  axval; 

end 

for  itheta  =  l:Ntheta 

spun (itheta, j2, figno)  =  axval; 

end 

jl  =  j2; 

end 

deltar  =  plotr(2)  -  plotr(l); 

deltas  =  spun (1, 3, figno)  -  spun ( 1 , 2  ,  figno)  ; 

axval  =  spun (2 , 2 , figno)  -  (plotr (1) *deltas/deltar) ; 

if (  axval  >  maxval  ) 
maxval  =  axval; 

end 

for  itheta  =  l:Ntheta 

spun (itheta, 1, figno)  =  axval; 

end 

end 

normalizer  =  1.0/maxval; 
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o\°  o\°  (D  CD  o\°  o\°  o\° 


for  figno  =  l:lastFigure 
for  jl  =  l:Nraxial 

for  itheta  =  l:Ntheta 

spun (itheta, j 1, figno)  =  normalizer*spun (itheta, j 1, figno) 

end 

end 

end 


case  B:  maxima 


lse 

zMax  =  1.0; 

if  (  (  Ntplot  ==  128  )  I  I  (  Ntplot  ==  256  )  ) 

verticalLabel  =  'maximum  rel.  intensity  (4x4  Avg . ) '; 
else 

verticalLabel  =  'maximum  rel.  intensity'; 

end 

for  figno  =  l:lastFigure 

ji  =  i; 

for  j2  =  2:Nraxial 
plotted  =  0.0; 

for  iplot  =  l:Ntplot 

value  =  intensity (j 1, iplot, figno) ; 

if (  value  >  plotted  ) 
plotted  =  value; 

end 

end 

for  itheta  =  lrNtheta 

spun (itheta, j2, figno)  =  plotted; 

end 

jl  =  j2; 

end 

deltar  =  plotr(2)  -  plotr(l); 

deltas  =  spun (1, 3, figno)  -  spun ( 1 , 2 , figno) ; 

axval  =  spun (2 , 2 , figno)  -  (plotr ( 1 ) *deltas/deltar) ; 

for  itheta  =  l:Ntheta 

spun (itheta,  1,  figno)  =  axval; 

end 

end 


color  map 
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map  =  zeros (25, 3); 


map (01, 1 ) 

= 

0.0000; 

map ( 02 , 1 ) 

= 

0.0400; 

map  ( 03 , 1 ) 

= 

0.0600; 

map  ( 04 , 1 ) 

= 

0.0800; 

map ( 05 , 1 ) 

= 

0.1000; 

map  (06, 1 ) 

= 

0.1200; 

map  ( 07 , 1 ) 

= 

0.1400; 

map (08, 1 ) 

= 

0.1600; 

map (09, 1 ) 

= 

0.1800; 

map  (10, 1 ) 

- 

0.2000; 

map (11, 1 ) 

= 

0.2018; 

map  ( 12 , 1 ) 

= 

0.2035; 

map (13, 1 ) 

= 

0.1936; 

map  (14, 1 ) 

= 

0.1837; 

map (15, 1 ) 

- 

0.1661; 

map (16, 1 ) 

= 

0.1485; 

map  (17, 1 ) 

= 

0.1232; 

map (18,1) 

= 

0.0979; 

map (19, 1 ) 

= 

0.0650; 

map (20, 1 ) 

= 

0.0320; 

map (21, 1 ) 

= 

0.0160; 

map (22 , 1 ) 

= 

0.0000; 

map (23, 1 ) 

= 

0.0000; 

map  (24, 1 ) 

= 

0.0000; 

map (25, 1 ) 

= 

0.0000; 

map (26, 1 ) 

= 

0.0000; 

map  (27, 1 ) 

= 

0.0000; 

map (28, 1 ) 

= 

0.0000; 

map (29, 1 ) 

= 

0.0000; 

map ( 30 , 1 ) 

- 

0.0000; 

map  (31, 1 ) 

= 

0.0000; 

map  ( 32 , 1 ) 

= 

0.0000; 

map ( 33 , 1 ) 

= 

0.0000; 

map  ( 34 , 1 ) 

= 

0.0000; 

map ( 35 , 1 ) 

- 

0.0000; 

map (36, 1 ) 

= 

0.0000; 

map  ( 37 , 1 ) 

= 

0.0000; 

map (38, 1 ) 

= 

0.0000; 

map (39, 1 ) 

= 

0.0000; 

map  (40, 1 ) 

= 

0.0000; 

map (41, 1 ) 

= 

0.0134; 

map  ( 42 , 1 ) 

= 

0.0269; 

map (43, 1 ) 

= 

0.1120; 

map  (44, 1 ) 

= 

0.1971; 

map (45, 1 ) 

= 

0.2899; 

map (46, 1 ) 

= 

0.3827; 

map  (47, 1 ) 

= 

0.4832; 

map (48, 1 ) 

= 

0.5837; 

map (49, 1 ) 

= 

0.6919; 

map ( 50 , 1 ) 

= 

0.8000; 

map (01,2) 

= 

0.0000; 

map (02,2) 

= 

0.0000; 

map  (03,2) 

= 

0.0000; 

map (04,2) 

= 

0.0000; 

map (05,2) 

= 

0.0000; 

map (06,2) 

= 

0.0000; 

map (07,2) 

= 

0.0000; 

map (08,2) 

= 

0.0000; 

map (09,2) 

= 

0.0000; 

map (10,2) 

= 

0.0000; 

map (11,2) 

= 

0.0000; 

map (12,2) 

= 

0.0000; 

map (13,2) 

= 

0.0000; 

map (14,2) 

= 

0.0000; 

map (15,2) 

= 

0.0000; 

map (16,2) 

= 

0.0000; 

map (17,2) 

= 

0.0000; 

map (18,2) 

= 

0.0000; 

map (19,2) 

= 

0.0000; 

map (20,2) 

= 

0.0000; 

map (21,2) 

= 

0.0296; 

map (22,2) 

= 

0.0493; 

map (23,2) 

= 

0.0976; 

map (24,2) 

= 

0.1459; 

map (25,2) 

= 

0.2019; 

map (26,2) 

= 

0.2579; 

map (27,2) 

= 

0.3216; 

map (28,2) 

= 

0.3853; 

map (29,2) 

= 

0.4567; 

map  (30,2) 

= 

0.5280; 

map (31,2) 

= 

0.5840; 

map (32,2) 

= 

0.6400; 

map (33,2) 

= 

0.6600; 

map (34,2) 

= 

0.6800; 

map (35,2) 

= 

0.7000; 

map (36,2) 

= 

0.7200; 

map (37,2) 

= 

0.7400; 

map (38,2) 

= 

0.7600; 

map (39,2) 

= 

0.7800; 

map (40,2) 

= 

0.8000; 

map (41,2) 

= 

0.8200; 

map (42,2) 

= 

0.8400; 

map (43,2) 

= 

0.8600; 

map (44,2) 

= 

0.8800; 

map (45,2) 

= 

0.9000; 

map (46,2) 

= 

0.9200; 

map (47,2) 

= 

0.9400; 

map (48,2) 

= 

0.9600; 

map (49,2) 

= 

0.9800; 

map  (50,2) 

= 

1.0000; 

map (01,3) 

= 

0.0000; 

map ( 02 , 3 ) 

= 

0.0077; 

map (03,3) 

= 

0.0192; 

map ( 04 , 3 ) 

= 

0.0307; 

map (05,3) 

= 

0.0499; 

map (06,3) 

= 

0.0691; 

map ( 07 , 3 ) 

= 

0.0960; 

map (08,3) 

- 

0.1229; 

map (09, 3 ) 

= 

0.1575; 

map (10,3) 

= 

0.1920; 

map (11,3) 

= 

0.2160; 

map ( 12 , 3 ) 

= 

0.2400; 

map (13,3) 

= 

0.2600; 

map (14, 3 ) 

= 

0.2800; 

map (15,3) 

- 

0.3000; 

map (16,3) 

- 

0.3200; 

map (17, 3 ) 

= 

0.3400; 

map (18,3) 

- 

0.3600; 

map (19, 3 ) 

= 

0.3800; 

map (20,3) 

= 

0.4000; 

map (21,3) 

= 

0.4200; 

map (22, 3) 

= 

0.4400; 

map (23,3) 

= 

0.4600; 

map (24, 3 ) 

= 

0.4800; 

map (25,3) 

= 

0.5000; 

map (26,3) 

= 

0.5200; 

map (27, 3 ) 

= 

0.5400; 

map (28, 3 ) 

- 

0.5600; 

map (29, 3 ) 

= 

0.5800; 

map (30,3) 

= 

0.6000; 

map (31,3) 

= 

0.5970; 

map (32,3) 

= 

0.5939; 

map (33,3) 

= 

0.5472; 

map ( 34 , 3 ) 

= 

0.5005; 

map (35,3) 

= 

0.4461; 

map (36,3) 

= 

0.3917; 

map (35,3) 

= 

0.3296; 

map (36,3) 

= 

0.2675; 

map (38, 3 ) 

= 

0.1978; 

map (40,3) 

= 

0.1280; 

map (41,3) 

= 

0.0640; 

map (42,3) 

= 

0.0000; 

map (43,3) 

= 

0.0000; 

map (44, 3 ) 

= 

0.0000; 

map (45, 3 ) 

= 

0.0000; 

map (46,3) 

= 

0.0000; 

map (47,3) 

= 

0.0000; 

map (48,3) 

- 

0.0000; 

map (49, 3 ) 

= 

0.0000; 

map ( 50 , 3 ) 

= 

0.0000; 

o_ 


183 


output  initial  pulse 


o, _ 

o 

disp(  'plotting...'  ); 

propagDist  =  0.0; 

figno  =  1; 

handle  =  figure (  1  ) ; 

pulseEnergy  =  energies (figno) ; 

etext  =  sprintf (  'pulse  energy  (FWHM)  =  %8.6f  J',  pulseEnergy  ); 

text  =  sprintf (  'prop,  distance  =0.0,  %s',  etext  ); 

axis (  [-radialMax  radialMax  -radialMax  radialMax  0.0  zMax]  ); 

caxis (  [  0.0  zMax]  ); 

title (  text,  ' FontWeight ' ,  'bold'  ); 

xlabel (  radialLabel,  'FontWeight',  'bold'  ); 
ylabel (  radialLabel,  'FontWeight',  'bold'  ); 
zlabel (  verticalLabel ,  'FontWeight',  'bold'  ); 

if (  get (  handles . draw3D,  'Value'  )  >  0  ) 
view (  [45  70]  )  ; 

else 

view (  [9090]  ) ; 

end 

hold  on; 

surf (  xx,  yy,  spun figno)  ); 

shading  interp; 

colormap (  map  ) ; 
colorbar ; 

if (  get (  handles . draw3D,  'Value'  )  >  0  ) 
grid  on; 

end 

hold  of f ; 

refresh (  1  ) ; 

if (  animation  ~=  0  ) 

frames (figno)  =  getframe (  handle  ); 

end 

O, _ 

o 

%  propagation  loop 

o, _ 

o 
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for  figno  =  2:Numplots 

propagDist  =  distances (figno) ; 

text  =  sprintf (  '  frame  %d  of  %d',  figno,  Numplots  ); 

set (  handles . status ,  'String',  text  ); 

if (  propagDist  >  0  ) 

pulseEnergy  =  energies (figno) ; 

if (  animation  ==  0  ) 

handle  =  figure (  figno  ) ; 

else 

elf  ; 

end 

if (  propagDist  <  (10A-6)  ) 

dtext  =  sprintf (  'prop,  distance  =  %5.3f  nm',  100000000 . 0*propagDist  ); 
elseif (  propagDist  <  (10A-3)  ) 

dtext  =  sprintf  (  'prop,  distance  =  %5.3f  micron',  1000000 . 0*propagDist 

)  ; 

elseif (  propagDist  <1.0  ) 

dtext  =  sprintf (  'prop,  distance  =  %5.3f  mm',  1000 . 0*propagDist  ); 
else 

dtext  =  sprintf (  'prop,  distance  =  %5.3f  m',  propagDist  ); 

end 

etext  =  sprintf (  'pulse  energy  (FWHM)  =  %8.6f  J',  pulseEnergy  ); 
text  =  sprintf (  '%s,  %s',  dtext,  etext  ); 

axis (  [-radialMax  radialMax  -radialMax  radialMax  0.0  zMax]  ); 

caxis (  [  0.0  zMax]  ); 

title (  text,  ' FontWeight ' ,  'bold'  ); 

xlabel (  radialLabel,  'FontWeight',  'bold'  ); 
ylabel (  radialLabel,  'FontWeight',  'bold'  ); 
zlabel (  verticalLabel ,  'FontWeight',  'bold'  ); 

if (  get (  handles . draw3D,  'Value'  )  >  0  ) 
view (  [45  70]  ) ; 

else 

view (  [90  90]  )  ; 

end 

hold  on; 

surf (  xx,  yy,  spun (:,:, figno)  ); 
shading  interp; 
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colormap (  map  )  ; 
colorbar ; 

if (  get (  handles . draw3D,  'Value'  )  >  0  ) 
grid  on; 

end 

hold  off; 

if (  animation  ==  0  ) 
refresh (  figno  ); 

else 

drawnow; 
refresh (  1  ) ; 

frames (figno)  =  getframe (  handle  ); 

end 

end 

end 

o, _ 

o 

%  save  movie 

o, _ 

o 

if (  animation  ~=  0  ) 

disp(  'saving  to  avi . . . '  ); 

text  =  sprintf (  'saving  to  AVI  (%d  frames) ',  Numplots  ); 
set (  handles . status ,  'String',  text  ); 
if (  animation  ==  2  ) 

movie2avi (  frames,  'spot.avi',  ' FPS ' ,  3,  'COMPRESSION',  'None'  ); 
else 

movie2avi (  frames,  'spot.avi',  'FPS',  3,  'COMPRESSION',  ' Indeo5 ' , 
'QUALITY',  100  ) ; 
end 

disp(  '  spot.avi  saved'  ); 

set (  handles . status ,  'String',  'spot.avi  saved'  ); 

end 

set (  handles .plot,  ' BackgroundColor ' ,  'green'  ); 

set (  handles .plot,  'Enable',  'on'  ); 

set (  handles . closeall ,  'BackgroundColor',  'green'  ); 

set (  handles . closeall ,  'Enable',  'on'  ); 

o _ 

O - 

%  respond  to  button  press  (close  all) 

o _ 

O - 

function  closeall^Callback (  hObject,  eventdata,  handles  ) 
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o\°  o\°  o\°  c/]  (D  (D  H-  o\°  o\°  o\° 


hObject  handle  to  closeall  (see  GCBO) 

eventdata  reserved  -  to  be  defined  in  a  future  version  of  MATLAB 
handles  structure  with  handles  and  user  data  (see  GUIDATA) 

f(  get (  handles . noanimation,  'Value'  )  >  0  ) 

Numplots  =  get (  handles . closeall ,  'UserData'  ); 

for  figno  =  l:Numplots 

handle  =  figure (  figno  ) ; 

close (  handle  ) ; 

end 

lse 

handle  =  figure (  1  )  ; 
close (  handle  ) ; 


et (  handles . closeall ,  'Enable',  'off'  ); 

END 
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